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LIST  OF  SYMBOLS 


y,z  Dimensionless  radial  and  axial  coordinates  (¥  =  ^ (y , z)  system) 

Coordinates  are  non-dimensionalized  with  respect  to  the  far  upstream 
hub  radius . 

Y,x  Dimensionless  radial  and  axial  coordinates  (Y  =  Y(x,¥)  system) 

R  Dimensionless  far  upstream  tip  radius 

G  Dimensionless  tip  radius  (G  -  G(x)) 

F  Dimensionless  hub  radius  (F  =  F(x)) 

¥  Dimensionless  streamf unction,  defined  by  Eqs  (2)  or  Eqs  (10) 

U,V,W  Dimensionless  velocity  components  in  radial,  azimuthal  and  axial 

directions.  Velocities  are  non-dimensionalized  with  respect  to  the 
far  upstream  axial  velocity. 

Dimensionless  vector  velocity,  components  (U,V,W) 

£  Dimensionless  angular  momentum,  £  =  YV  =  yV 

M  Mach  number 

p  Dimensionless  density.  Density  is  non-dimensionalized  with  respect  to 

the  far  upstream  static  density. 

T  Temperature 

S  Entropy 

Specific  heat  at  constant  volume 

H  Dimensionless  stagnation  enthalpy.  Enthalpies  are  non-dimensionalized 

with  respect  to  the  square  of  the  far  upstream  axial  velocity. 

h  Dimensionless  static  enthalpy 

P  Dimensionless  pressure.  Pressure  is  non-dimensionalized  with  respect 

to  the  product  of  the  far  upstream  static  density  and  square  of  the 
far  upstream  axial  velocity. 

K  Group  of  term1:  defined  by  Eq  (15) 

V  Variational  functional,  defined  by  Eq  (14) 

ft  Dimensionless  rotor  angular  velocity.  Non-dimensionalized  with  respect 

to  the  far  upstream  axial  velocity  divided  by  the  far  upstream  hub 
radius. 
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SUBSCRIPTS 

d  Refers  to  (jump)  conditions  at  a  disc 

o  Conditions  far  upstream 

e  Conditions  at  exit 

T  Stagnation  conditions 

~“,,°0  Conditions  far  upstream,  far  downstream 

NOTE:  Symbols  used  in  finite  element  analysis  and  numerical  analysis  are 

introduced  and  defined  within  the  text. 
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SECTION  I 


INTRODUCTION 


The  throughflow  theory  of  flow  in  turbomachines  describes  the  overall  ef¬ 
fect  of  an  entire  blade  row  or  rows  upon  the  fluid  properties  within  the  ma¬ 
chine.  The  throughflow  field  is  considered  to  be  axially  symmetric,  and  because 
of  its  position  as  the  "parent"  flow  field  to  the  cascade  field  and  the  second¬ 
ary  field,  it  is  important  that  the  throughflow  field  be  accurately  described. 

An  extensive  description  of  throughflow  theory  for  turbomachines  is  given 
in  Refs.  1-3.  In  Ref.  3,  a  variational  formulation  of  the  incompressible 
throughflow  problem  suitable  for  the  description  of  highly  rotating  flows 
existing  in  annuli  with  large  variations  in  hub  and  tip  radii  is  described. 

The  present  study,  described  herein,  extends  the  study  of  Ref.  3  to  include 
the  effects  of  compressibility.  The  effects  of  entropy  production  within  the 
blade  rows  is  included  in  the  formulation  though  the  flow  external  to  the  blade 
rows  must  be  considered  perfect.  In  addition  the  meridional  Mach  number  must 
be  less  than  unity. 

This  latter  restriction  is  not  very  limiting  with  regards  to  the  extent 
of  applicability  of  the  results,  because  few  (if  any)  turbines  or  compressors 
are  contemplated  that  would  utilize  meridional  Mach  numbers  in  excess  of  unity. 
Of  course,  the  full  flow  Mach  number  (including  the  swirl  component  of  velocity) 
may  be  much  in  excess  of  unity.  It  should  be  noted,  however,  that  it  is  the 
approach  to  unit  meridional  Mach  number  that  is  at  the  heart  of  the  computa¬ 
tional  difficulties  of  any  program  that  calculates  near-sonic  flow  fields. 

Any  iterative  program  utilizes,  somewhere  in  the  program,  a  "guess"  for  tne 
streamline  location  that  will  later  be  adjusted  to  numerically  satisfy  the 
describing  equations.  A  glance  at  subsonic  flow  tables  (Ref.  4  for  example) 
shows  that  a  one  percent  error  in  streamtube  area  would  lead  to  the  prediction 
of  sonic  flow  for  a  flow  whose  "true"  Mach  number  should  be  0.9.  Without 
special  corrective  measures  being  taken,  this  error  would  lead  to  numerical 
divergence  of  the  density  calculation. 

There  are  two  reasons  for  the  comparatively  great  difficulties  encountered 
in  the  numerical  calculation  of  highly  rotating  compressible  flow  fields.  The 
presence  of  the  large  rotation  leads  to  highly  distorted  axial  velocity  pro¬ 
files,  and  the  presence  of  the  large  swirl  component  of  velocity  leads  to  large 
total  flow  Mach  numbers.  The  fluid  density  is,  of  course,  determined  by  the 
total  flow  Mach  number  so  the  combination  of  axial  velocity  distortion  and 
large  density  variations  make  it  difficult  to  estimate  an  accurate  "first 
guess"  for  streamline  position.  As  a  result  erroneous  "supersonic"  portions 
of  the  flow  field  can  easily  appear  in  the  initial  portions  of  the  calculations, 
even  in  a  flow  field  that,  when  eventually  determined,  will  have  rather  modest 
meridional  Mach  numbers  throughout. 

These  preliminary  remarks  are  intended  to  alert  the  reader  to  the  impor¬ 
tance  of  wisely  selecting  the  iterative  procedures  of  calculations  such  as 
those  described  herein,  and  to  explain  why  such  an  emphasis  has  been  placed 
upon  the  method  and  order  of  the  several  subsidiary  calculations  involved  in 
the  program. 
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SECTION  II 


VARIATIONAL  FORMULATION  OF  THE 
COMPRESSIBLE  THROUGHFLOW  PROBLEM 


1.  FORMULATION  OF  THE  EQUATIONS 


The  pertinent  equations  are  developed  in  detail  in  Ref.  3,  but  for  con¬ 
venience  non-dimensional  forms  of  the  resultant  equations  will  be  repeated 
here.  When  perfect  flow  is  considered,  conservation  of  angular  momentum  and 
Euler's  momentum  equation  indicate  that  external  to  the  blade  rows  the  angular 
momentum  and  stagnation  enthalpy,  respectively,  are  conserved  along  stream- 
surfaces  . 

Information  regarding  variations  across  streamsurf aces  is  obtained  by  con¬ 
sidering  the  normal  component  of  the  combined  momentum  and  thermodynamic  equa¬ 
tions.  The  resultant  equation  gives  a  relationship  for  the  tangential  vorti- 
city  which  may  be  written  in  the  form: 


3_(  1_  34\  3,  1_ 

3z  py  3z  '  3y  ^  py  3y  ' 


3S  _  1  3£2i 

3<F  "  2y2  3V  J 


(1) 


This  form  is  valid  for  flow  external  to  the  blade  rows,  provided  only  that 
the  flow  external  to  the  blade  rows  may  be  considered  nonviscous. 

The  dimensionless  velocity  components  are  related  to  the  dimensionless 
streamf unction  by: 
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y  3y 


(2) 


The  boundary  conditions  on  Y  are 

1  R 2 

V  =  -  on  y  =  F(z)  ,  T  =  -  —  on  y  =  G(z)  (3) 

When  a^'  <ator  discs  are  considered  to  exist  in  the  flow  field,  the  matching 
conditions  across  the  discs  follow  from  the  continuity  equation  and  radial 
momentum  (assuming  no  radial  forces)  to  give  respectively  the  jump  conditions 
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In  general  there  will  be  jumps  in  tangential  momentum  and  stagnation  en¬ 
thalpy  imparted  by  the  blade  rows,  which  must  be  related  by  the  Euler  momentum 
equation  to  give 


[H]d  =  nmd  (6) 

Finally,  the  thermodynamic  equation  in  conjunction  with  the  equation  for 
the  (known)  stagnation  enthalpy  gives  the  equation  for  the  density. 


H  -  \  (U2  +  W2  +  = 

Utilizing  Eq  (2),  this  may  be 

(  2H  -  ^-)  p2  -  2h0pY+1  e 


rearranged  to  give: 
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Eqs  (1)  -  (8)  represent  the  mathematical  statement  of  the  problem.  Fig.  1 
indicates  the  various  geometric  relationships.  It  should  be  noted  that  any 
entropy  variations  occurring  are  considered  to  have  been  created  within  the 
blade  rows,  because  Eq  (1)  has  not  included  the  effects  of  viscous  terms. 

This  being  the  case,  external  to  the  blade  rows  the  entropy  is  a  function  of 
streamfunction  only.  This  model  gives  a  good  approximation  to  the  throughflow 
flow  field, (1)  and  in  this  report  only  the  special  case  of  perfect  flow  (S  = 

SQ  throughout)  is  actually  calculated.  The  formulation,  however,  is  not  re¬ 
stricted  to  isentropic  flow  fields. 


2.  FORMULATION  OF  THE  TRANSFORMED  EQUATIONS 


The  numerical  solutions  to  follow  are  conveniently  obtained  by  formulating 
the  equations  to  solve  for  the  radial  position  as  a  function  of  axial  position 
and  streamfunction.  We  thus  transform  the  independent  variables  from  y  and  z 
to  ¥  and  x,  where  x  =  z.  We  also  write  y  =  YCx.'F)  in  order  to  facilitate 
simple  conceptual  separation  of  the  two  systems.  Routine  transformation  of  the 
appropriate  equations  leads  to  the  following  set  of  equations  in  the  x , T  system. 

Differential  equation 
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Velocity  Relationships 
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Boundary  Conditions  (enclosed  flows  only) 
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Density  Relationship 


(  2H  -  £ )  p2  -  2h0pYJ'1  e  Cy 


(y||)2 


l  I 


The  Euler  equation  remains  as  in  Eq.  (6). 


3.  VARIATIONAL  STATEMENT  OF  THE  PROBLEM 


We  now  define  a  variational  functional,  V,  by 


V  *  .1  ;f  tp  ~  pe  +  P^2  +  W2}jydydz 


This,  of  course,  is  simply  the  integral  of  the  meridional  momentum  over 
all  space.  For  convenience  in  relating  the  compressible  results  to  the  in¬ 
compressible  results  we  define  the  parameter  K  by: 

K  "  p  "  pe  +  \  P  U2  0 


The  thermodynamic  equation  may  be  written 


U2  1 

TVS  -  VH  -  V  %  -  -  VP 
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so  that  with  Eq  (15), 

TT2 

VK  -  p{VH  -  TVS}  +  ~~  Vp  (17) 

Utilizing  the  definitions  of  Eqs  (14)  and  (15),  we  obtain  the  expression 
for  V  in  the  (x,^)  system  as: 

oo  -R2  /2 

?  -  /  /,  [  K  +  |  {  U2  +  W2  -  V2}]  Y  ~  d'i'dx  (18) 

-00  -'/2 

The  variation  of  this  quantity  is  now  taken.  Noting  from  Eq  (17)  that 


SK  -  p{<5H  -  T6s}  +  —  6p  (19) 


the  variation  of  Eq  (18)  utilizing  Eqs  (10)  leads,  after  some  manipulation, 
to 


6V 


00  /  2  JNy  <“.y  1 

/  /  [  p{6H  T6S}  Y  +  K  6(Y  -^)  +  ~  6 

’  -V2 


1  +  (f ) 

9x 


9Y 

9Y 


"  2  6{  Y~  H  dVFdx  (20) 


In  the  flows  considered  in  this  report,  perfect  flow  Is  assumed  to  exist 
outside  the  blade  rows.  This  being  the  case,  both  the  entropy  and  stagnation 
enthalpy  are  constant  along  streamsurfaces .  When  the  variation  is  taken  in  the 
Xj'J'  system,  as  it  is  here,  the  variation  of  functions  of  ¥  alone,  is  of  course, 
zero.  It  follows  that  neither  5H  nor  6S  contribute  in  Eq  (20).  Utilizing  this 
fact  and  rearranging  the  terms  in  Eq  (20)  we  then  obtain: 


«  -R2/2 

6v  =  ;  / 

-00  -1/2 


6Y 


l?([KY 


2  Y 


dK 

i  9pii 

_  9_ 

9Y 

9x  Y 

9 

1  + 

<-)2 

9Y 

2Y  9¥ 

9x 

PY  —  2 

Pi  9Y 

9Y 

P 

(YiI)2 
'3  r 

1  1 

♦ 

2 

_  Y 

6y]  + 

9x 

9Y 

9x 

5Y  " 

d¥dx 

1 

2 

P 

PY 

9Y 

9¥  ' 

(21) 


Eq  (17)  with  Eqs  (10),  gives  us 


9K  ,  ,3H  cS 

9T  "  p  l3¥  9Y 


1 

2 


(-) 


9Y  2 
(Y  —  ) 


il_  la 

2Y2  9Y 


(22) 
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J 


Substituting  this  expression  into 
terms  leads  to 


Eq  (21), 


and  integrating  the  last  two 


+ 


dY  (23) 

+oo 


Interior  to  the  domain,  6Y  is  arbitrary.  It  can  be  seen  from  Eq  (9)  that 
the  first  term  gives  the  Euler-Lagrange  equation  governing  Y(x,'t').  Utilizing 
Eqs  (10)  and  (15),  it  is  apparent  that  the  second  term  vanishes  if  the  boundary 
is  specified  (6Y=0  on  =  -1/2,  R2/2)  or  if  the  boundary  is  a  free  streamline 
(P=Pe),  The  remaining  terms  arise  from  the  far  upstream  conditions,  jump  con¬ 
ditions  at  the  actuator  discs,  and  far  downstream  conditions.  For  enclosed 
flows,  radial  equilibrium  at  far  upstream  and  far  downstream  ensures  that 
9Y 
“  gx 

U  =  - ^  is  zero.  Eqs  (12)  indicate  that  the  matching  conditions  at  the 

pY  av 

actuator  discs  are  also  naturally  satisfied  by  this  variational  statement. 

It  can  now  be  seen  that  the  variational  principle  represents  a  complete 
statement  of  the  physical  problem  in  which  the  effects  of  compressibility  and 
entropy  variations  have  both  been  included.  The  boundary  conditions  and  all 
matching  conditions  are  all  natural  conditions  of  the  variational  formulation. 
As  a  result,  a  finite  element  variational  approach  may  be  employed  without  pre¬ 
scription  of  these  boundary  and  matching  conditions,  and  they  will  auto¬ 
matically  be  satisfied  within  the  numerical  approximation. 
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4.  FORMULATION  OF  THE  VARIATIONAL  EQUATION 


It  still  remains  to  represent  6V  in  terns  of  prescribed  variables  (H,£) 
and  functions  of  Y  and  6Y.  Eq  (20)  and  the  discussion  immediately  following 
Eq  (20)  lead  to 


6Y 


‘  2  6  {Y~W}]  d4/dx 


(24) 


We  now  restrict  our  attention  to  the  calculation  of  isentropic  flow  fields. 
An  expression  for  K  is  obtained  by  first  noting  that  the  dimensionless  pressure 
may  be  written 


PYP, 


Utilizing  Eq  (13)  with  S  =  SQ,  it  follows  that 


(25) 


(26) 


Thus  with  Eqs  (10),  (15)  and  (26)  it  follows  that 


1  +  1 

1  3Y 
1  3x 

)* 

Y  ■ 

3Y1 

34'J 

2 

+  P 


li" 

y2_ 


- 


(27) 


The  desired  form  of  the  functional  variation  is  then  obtained  to  give: 


(28) 


8 


The  subsidiary  equation  for  the  density,  for  the  case  of  isentropic  condi¬ 
tions  is: 


|»  -  £]  -  V" 


1  + 

3Y 

3x 

f 

Y  III 

1  3fJ 

(29) 


Usually  the  far  upstream  approach  Mach  number  will  be  prescribed,  so  we 


note 


H0  -  h0  -  j  -  hQ  ( ~  M02) 


hence  h„ 


JL _ 1 

Y-l  M  2 


(30) 


and  H„ 


i  1  +  2rMo! 


Y-l 


V 


(31) 


The  local  stagnation  enthalpy  H  follows  from  application  of  the  Euler 
Momentum  equation,  Eq  (6). 

It  is  important  to  note  that  the  density  variations  were  included  in  the 
variation  of  functional  V  in  Eq  (18)  to  determine  6V  in  Eq  (20).  However, 
collecting  terms  involving  6p  produces  the  governing  differential  equation 
multiplying  6p  as  the  density  variation  contribution.  Then,  since  the  govern¬ 
ing  equation  is  satisfied  as  the  Euler  equation  from  consideration  of  varia¬ 
tions  6Y,  the  expression  involving  6p  is  zero.  That  Is,  we  may  write 

*  S  5Y  +  f  6P  "=  I'  6Y  "  0 


and  the  variational  problem  reduces  to  satisfying  Tjy  6Y  =  0.  Then  the  finite 

element  approximation  need  introduce  only  an  approximation  to  Y(x,y)  in  terms 

i  3  V 

of  solution  point  values  Y.  and  the  variational  problem  requires  — r  =  0  sub- 

J  9YJ 

ject  to  the  prescribed  boundary  conditions. 

The  unknown  density  field  still  appears  in  the  coefficients  of  the  varia¬ 
tions  of  position  terms  and  satisfies  the  subsidiary  equation,  Eq  (29).  Thus 
given  any  streamfield,  the  associated  density  field  is  completely  specified. 

These  two  calculation  steps,  (a)  the  streamfield  position,  and  (b)  the 
corresponding  density  field, constitute  the  central  feature  of  the  coupled 
iterative  scheme  described  in  the  following  sections. 
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We  remark  that,  as  an  alternative  scheme,  the  density  field  could  be  in¬ 
cluded  as  an  unknown  in  the  variational  problem  itself,  by  expanding  '{(*,'¥)  = 

fj(Yj)  and  PCx.f)  *  f2(p*).  The  formulation  and  computation  become  complicated, 

since  the  subsidiary  density  equation  must  be  imbedded  as  a  constraint  by 
Lagrange  multiplier  methods.  The  vector  of  finite  element  nodal  unknowns  is 

thereby  increased  to  include  (Y*,  pj,  X*},  where  X(x.’l')  is  the  Lagrange  multi¬ 
plier.  In  turn,  this  increases  the  vector  of  unknowns  in  the  nonlinear  alge¬ 
braic  system  resulting,  and  the  problem  formulation  is  computationally  inferior 
to  that  of  the  coupled  position  and  density  problems  above. 


5.  INTERPRETATION  OF  TERM  K 


It  is  clear  that  the  effect  of  finite  Mach  number  is  felt  directly  in  the 
coefficients  of  the  last  two  terms  of  Eq  (28)  through  its  effect  on  the  dimen¬ 
sionless  density  p.  It  is  of  interest  to  consider  also  the  effect  of  Mach 

Y-l  Pt 

number  on  the  term  K.  Writing  -rr-  H  =  — ,  we  have  from  Eqs  (10)  and  (27) 

Y  Pt 


K 


or  K 


r.  2 


p  -p.-e^n 


2Y 


Ii  1  *1*’ 

°T  1  +  ^  M* 


-  P 


1  +  |  M2 

- - - _L 

(1  +  ~  M2}  Y-l 


-  P. 


We  note  upon  expansion  that  K  = 


P  -  ~  M4P  + 
e  8  T 


Thus  for  small  Mach  numbers,  the  compressibility  effect  upon  the  parameter  K 
is  very  small  indeed. 


SECTION  III 


FIN' YE  ELEMENT  METHOD 


1.  INTRODUCTION 


The  preceding  analysis  develops  a  variational  formulation  of  the  through- 
flow  problem.  The  Euler  equation  is  the  governing  partial  differential  equa¬ 
tion  for  the  throughflow,  and  natural  boundary  conditions  prescribe  conditions 
at  the  actuator  discs  and  at  downstream  infinity.  For  the  present,  only  en¬ 
closed  flows  are  considered.  In  this  analysis  the  f ully-inf inite  flow 
domain  is  idealized  as  a  finite  but  extensive  region  with  boundary  conditions 
at  infinity  applying  on  the  remote  boundaries  x  =  x0  and  x  =  xjg.  The  finite 
domain  approximation  is  not  essential  to  the  method.  Since  the  asymptotic 
far  field  solution  is  known,  this  may  be  inclined  in  the  finite  elenent  analysis 
as  an  "exterior"  element  approximant  matching  the  solution  at  the  remote  bound¬ 
aries  x q  and  Xfl.(s)  However,  it  is  simpler  aid  adequate  to  use  the  finite 
domain  approximation  and  the  natural  boundary  condition  insures  parallel  flow 
as  the  solution  approaches  the  radial  equilibrium  condition  at  the  ends. 

It  was  remarked  in  the  variational  development  that  the  inverse  formulation 
in  Y (x,f)  achieved  by  transformation  to  (x,’!')  coordinates  is  particularly  ad¬ 
vantageous.  There,  it  was  pointed  out  that  enthalpy  II  and  angular  momentum  £ 
are  constant  along  streamsurfaces  in  the  flow  between  discs.  In  addition,  the 
domain  is  mapped  to  a  rectangular  region  bounded  by  T  =  ,  =  -1/2  and 

R2 

¥  =  'Fj- £p  =  -  —  streamlines  and  with  intermediate  ¥  lines  parallel  to  the  bound¬ 
ing  streamlines.  A  representation  of  rectangular  elements  is  a  natural  choice 
in  the  (x,4')  plane.  There  is  no  domain  error  at  the  hub  and  tip  and  the  use  of 
parametric  element  forms,  necessary  in  an  accurate  (x,y)  plane  approximation  at 
the  walls,  is  thereby  avoided. 

The  finite  element  description  of  the  problem  is  now  developed  for  a 
general  element  form.  Following  this,  detailed  derivation  is  presented  for  a 
low  order  element  approximant,  the  bilinear  interpolant. 

2.  NATURE  OF  THE  FINITE  ELEMENT  SYSTEM 


The  desired  variational  functional,  V,  was  described  in  the  previous 

3  f  3  f 

section.  We  now  utilize  the  notation  =  fy,  =  fx,  f  any  function  to  re¬ 
write  the  variation  6Y  as  given  in  Eq  (28)  as: 

'  1  J2  [[^PH  +  ^tpI^+0$)]s<Y^>+ssli^7') 


<5Y  = 


'  2  6  (  Y~  YV 


d^dx 


(32) 
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It  is  to  be  noted  that  the  term  Pg  has  been  excluded  for  algebraic  simpli¬ 
city.  In  the  situation  considered  here  of  enclosed  flows  only,  Pe  makes  no 
contribution  to  the  flow  behavior. 

Assume  Lagrangian  interpolation  on  any  given  mesh  of  elements  and  let  the 
nodes  be  numbered  consecutively.  For  example,  on  a  rectangular  mesh  of  size 
S  ■  (MxN)  the  general  node  k  defined  by  the  ordinate  index  pair  (i,j)  is  k  = 
(i-l)N  +  j.  The  finite  element  approximation  may  be  written  in  terms  of  the 

nodal  values  {Y  },  k=l,  ...fS,  as 

Y(x,t)  =  Yk4>k(x,4')  03) 


where  the  repeated  index  implies  summation. 

The  functions  ^(x*^)  are  termed  nodal  or  global  basis  functions.  They  are 
non-zero  only  over  those  elements  immediately  adjacent  to  node  k,  assume  unit 
value  at  node  k  and  are  zero  at  all  other  nodes. 

Substituting  for  Y(x,'f)  into  V(Y),  then  V(Y)  becomes  V (Y  ,  { } ) ,  with 


6Y 


-  \ 
oi 


0 


which  implies 


— r-  =  0  ,  k=l ,  .  .  .  ,S 

3Y 


(34) 


This  represents  the  system  of  finite  element  equations  for  the  discretized 
throughflow  problem.  From  Eq  (32),  the  detailed  expression  for  Eq  (34)  is 


(35) 
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We  examine  the  expression  — r  (YYy)  in  detail: 

ox 


<Y°VSy 


A  ^Vw 

oY 


6.  Ye<f>  (4>q)  +  6.  QYa<J>  (<J>0) 

ka  TaVYB  \j/  kg  Ya  8  y 

Y%,  (4>0)  +  Ya4>  (<J).  ) 

Tk  Yp  f  Ya  k  4/ 


or 


“IT  <YV 

3Y 


Ya{(|),  (0  )  +  $(♦,)' 

’k  Yf  Ya  k  4/ ■ 


-J_| 

’  3Ykl 


|l  +  Yx2 

and  — 

[S2  YT» 

l  YY 

and  . 

Y  ] 

The  other  terms  of  interest 


may  be  similarly  expanded.  Substitution  into  Eq  (35)  and  integration  over  the 
"patch"  of  elements  connected  at  node  k  determines  the  general  finite  element 
equation  corresponding  to  node  k. 


The  additional  variables,  enthalpy  H (¥)  and  angular  momentum  £(¥),  are  known 
functions  throughout  the  flow  field;  the  density  P  is  an  implicit  function  of 
the  streamline  field  and  is  assumed  known  corresponding  to  the  prior  position 
iteration.  Tills  completes  the  most  general  finite  element  description  of  the 
problem. 


3.  THE  ELEMENT  APPROXIMATIONS 

In  the  above  statement,  a  very  general  nodal  or  global  basis  {4>^(x,4')}  was 
introduced.  The  problem  reduces  to  quadrature  on  each  nodal  patch.  In  turn, 
this  quadrature  is  most  conveniently  evaluated  as  an  accumulated  quadrature 
from  each  element  of  a  given  patch.  To  do  so,  the  global  basis  {^(x,'F)}  is 
resolved  to  a  local  element  basis  on  each  element. 

Consider  a  partition  of  the  flow  .  ield  by  rectangular  elements.  The 
simplest  approximation  uses  bilinear  interpolation  on  each  element  and  is  de¬ 
veloped  here.  The  extension  to  higher  order  elements  is  immediate.  A  general 
element  and  nodal  patch  with  global  basis  function  is  shown  in  Fig,  2. 

For  convenience  of  notation  in  deriving  the  finite  element  equations, 
replace  node  k  by  the  ordinate  node  pair  (m,n)  as  market  with  k  -  (m-l)N  +  n. 

Let  x  .  -  x  =  a  and  ¥  .  -  ¥  =  b  and  introduce  the  normalized  element 

nt+1  m  n+1  n 
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Element  Approximants 


coordinates  x  =  (x  -  x  )/a  and  ?  =  (T  -  4*  )/b.  In  bilinear  interpolation  on  an 

m  n 

e’eraent  e,  Y  (x.'F)  =  a  +  a  ,  x  a,  4*  +  a.xY  with  the  {a,}  determined  such  that 
e  o  i  2  3  i 

Y  (x.'F)  interpolates  the  four  corner  values  Ym,  Ym+''',  Ym, ,  ,  and  Y11^.  In 
e  n  n  n+1  n+1 

normalizeu  coordinates,  the  bilinear  inter)  olant  is 

Y  (x ,?)  =  (1  -  x)(l  -  40  Y1"  +  5(1  -  ?)  Y1^1  +  xY  y"*1  +  (1  -  x)?  Y®  (36) 
e  n  n  n+1  n+l 

The  statement  of  Eq  (34)  may  be  rephrased  in  terms  of  element  contributions. 

The  functional  Y  is  represented  as  a  sum  of  element  quadratures  by 


M-l  N-l 

V  =  E(Y)  =  Z  E  Ym 
e  e  _ ,  n 


where  Y  =  Ym  is  the  contribution  to  Y  of  element  e  identified  by  the  lower- 
e  n 

left  index  pair  (r,n). 


Eq 


(34)  becomes  for  any  node 


3  Y  9  Ye 

3y!  e  3Y7- 
J  j 


M-l  N-l 

z  z 

m=  1  n=  1 


(i,j), 


n 


(37) 


Consider  the  contribution  of  a  single  element  at  (m,n)  to  the  system.  Then, 
writing, 


K(x,4) 


Y-l 

Y 


PH 


+ 


the  element  contribution  is 


9  Y 
3Y 


r  nr+ 

■1  r  n+1  - 

J 

X 

V 

K(x,Y) 


3Y1 

J 


(YY^) 


+ 

2p  3Y1 
J 


1  +  Y" 


YY„ 


-)  -  - 


P  _jL 

2  3Y7 
J 


(f-  Yy)  ]  dYdx 


(38) 


In  the  interests  of  developing  as  simple  a  numerical  problem  as  possible 
for  this  initial  analysis,  and  motivated  by  the  analogous  assumption  that 

=  constant  for  the  incompressible  calculations ,( 3 )  the  quadrature  problem  is 

simplified  by  the  assumption  K(x,¥)  =  =  constant  in  an  element.  This 

approximation  is  involved  only  in  the  quadrature.  The  enthalpy  and 
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angular  momentum  are  constant  along  f  lines  in  a  given  element.  For  the  quad¬ 
rature  He  =  (H^  +  ^n+j)/2  and  =  (£*  +  £*+^)/2  are  assumed  constant  in  an 


element.  Similarly,  the  density  p  is  constant  at  the  mean  value 
=  l/4{  P  m  +  p 


nH-1  m  mfl-, 

n  +  P„+l+ 


e  n 

Under  these  assumptions  the  quadrature  of  Eq  (38)  becomes 


m 


3V 
_ n 

3Y1 

J 


. 

e  3y: 


x  + 

mfl  n+1 


x 

J  m 


/  (YY  )  dYdx 

If  * 


till 

2  P  3Y1 


Xmfl  ^n+1  ,  1  +  Y2 

/  / 


¥ 


\  YY 


J  ™ 


¥ 


2  c2  X 

-jd^dx  ~  P  f 


Y 

mfl  n+1 

/  *  Y 

x  Y  ~Y  dH,dx 

m  n 


(39) 


The  three  integrals  in  Eq  (39)  are  evaluated  in  turn.  Consider  the  first 
integral , 


X  .  .  Y  4.1 
mfl  n+1 


I,  =  /  / 

x  Y 
m  n 


(YYy)  dWdx 


x 

Vl  n+1 


/  /  (  \  Y2 )y  dYdx 


x  Y 
m  n 


.  mfl 

-j  /  (Y2(x,Yn+1)  -  Y2(x,^)}  dx 

x 

m 


Applying  the  trapezoidal  rule, 


~  I  [(Ym  V 

4  U  n+l' 


(Y^1)2 

v  n+r 


-  - 


(Vt^1)2](. 


'mfl 


m 


(40) 


In  the  second  integral,  since  the  element  interpolant  is  bilinear,  the  ap- 
proximant  varies  linearly  along  the  element  sides.  Derivatives  Y^  and  Y^  at 

the  corner  nodes  can  be  replaced  by  the  corresponding  difference  quotients  in 
the  4-point  quadrature  formula 
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■J 


1C  ¥ 

nrfl  n+1 

f  /  f(x,¥)  d¥dx 


m 


¥ 


-  j  [fm  +  f11^1  +  f“  +  f“t}](¥  ..-¥  )(x..-x  ) 
Ann  n+1  n+1  n+1  n  m+1  m 


to  give 


n  + 


r2  si 


mfl  m 
n+1  ~  n+1 


x  -  x 
m+1  m 


m+1  nrfl  _  nrfl 
-  n+l\  n+1  n 


1  + 


nrfl  _  m 
n+1  n+1 


x  , -  x  . 
m+1  m  / 


1  + 


fm+l  _  Y«n 


x  .  1  -  x 
nrfl  m 


Ym  I 

n+1  n+1  “  n 


ym+l  ym+l  _  ym+lj 
n  n+1  n  I 


1  + 


ynrfl  -  Y 
n  n 


x  , ,  -  x 
nrfl  n 


vm /y®  vmN 

nk  n+1  n; 


(V,^n  )2(X..-X  )  (41) 
n+1  n  m+1  m 


Similarly,  the  third  quadrature  gives 


I.  ~  T 


rQ2  -  02 

nrfl  m+1 
n+1  n 


(ym  2_  (ym  2 
n+1  n 

Ym  ym 

n+1  n 


(x. 


nrfl 


O 


(42) 


Writing  I  =  I  (Ym,  Y®  ,  y"*1,  Y^J")  ,  a  =  1,2,3  ,  then  the  following 
a  a  n  n+1  n  n+1  & 

permutation  relations  are  evident  on  inspection 


T  fv®  Ym  Vmt  ^  — 

a  '  n  ’  n+1  ’  n  *  n+r 


a  Un+1  *  n  *  n+1  ’  n  ; 


=  I  (Y^1  ,  Y®^  ,  Ym  ,  Y®)  =  -  I  (y”*1  ,  Y"*1  ,  Ym 
a  n  n+1  '  n  n+1  a  n+1  n 


Then  differentiating  1^  with  respect  to  Y  we  have, 


n+1  * 


Ym) 

n 


(43) 
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The  term  -  '-.an  be  evaluated  and  - 

3Yin  3Ym 

n  n+ 


3Ym  ’  Sy"*1  ’  Sy"^ 
n+1  n  n+1 


all  follow  directly 


from  the  permutation  relations  in  Eq  (43). 

3Ia 

The  integrated  expressions  — r  ,  a  =  1,2,3  may  be  substituted  into 

3y! 

m  .1 

3Ym 

— r  .  We  write 


-in  m  m  „m+l  „m+l 1 
n  n  ’  n+1  ’  n  *  n+1, 


££i  + 

1 

31 2 

3Ym 
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3Ym 

r2 

6  3Ym 

Then, using  the  permutation  relations 


pm  Ym  y11^  ^  yinl  ^ 

n  n  ’  n+1  ’  n  ’  n+1 


"  F*( 
n| 


m  m  m+1 

n+1  ’  n  ’  n+1 


■C1) 


+  Fm[YnT+1  ,  Y"t;  ,  Ym  ,  Y"  ,]  -  f"/y”1  ,  Y1-1  ,  Y"  Yn|  (45) 

n\  n  n+1  n  n+lj  nj  n+1  n  n+1  n  / 

Eq  (45)  is  the  contribution  of  the  general  element  e  to  the  finite  element 
equation  at  node  (i,j).  The  product  of  delta  functions  in  F™  ensures  that  the 

correct  non-zero  element  contributions  are  accumulated.  For  example,  if 
(i,j)  =  (m+l,n)  then  the  contribution  of  the  element  is 


Fm(Ynrfl  ,  Y^J  ,  Ym  ,  Ym 
nl  n  n+1  n  n+1 


Individual  element  contributions  are  now  collected  to  determine  the  finite 
element  system. 


3  V 

ll  — " 
m  n  3Y1 
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An  arbitrary  internal  node  point  (m,n)  has  four  adjacent  elements  and  col¬ 
lecting  element  contributions  at  this  node, 


0  = 


3Ym 
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|  +  pm  ym  ym  y^1  Ym+lj  ( 
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At  a  nodal  point  (m,n)  on  the  boundary  only  2  elements  are  adjacent  except  at 
the  region  corners  where  there  is  a  single  element.  For  example,  at  a  point 
(m,n)  on  the  lower  boundary, 


,m-l  ym  +  m  Ynr+1  Ym+1  ym  m' 

1  *  2  *  1  |  112  ’  1  ’  2  ’  1 


0 


A  similar  relation  holds  on  the  upper  boundary  and  ai  the  remote  upstream  and 
downstream  boundaries. 

The  boundary  conditions  are  now  included  in  the  formulation  of  the  dis- 

i  N 

crete  problem.  The  enclosed  flows,for  instance, imply  that  Yj  and  Y,  are 

prescribed  for  all  i  and  j,and  these  values  are  set  directly  into  the  nonlinear 
system.  The  remote  end  conditions  are  retained  as  the  natural  boundary  condi¬ 
tion  resulting  from  arbitrary  end  node  variations,  and  satisfy 
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(47) 


Eqs  (46),  (47)  and  the  prescribed  boundary  conditions  (known  position)  determine 
the  nonlinear  system  of  finite  element  equations  for  streamline  position  for  a 
prescribed  density  field.  The  nonlinear  system  may  be  solved  and  the  new  stream- 
field  used  to  improve  the  density  field  solution.  The  coupled  solution  of 
streamline  position  and  density  is  then  iterated  to  convergence. 
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SECTION  IV 


NUMERICAL  ANALYSIS 

1.  FORMULATION  OF  COUPLED  ITERATIONS 


The  streamline  position  Y(x,4')  and  the  density  p(x,V)  are  the  primary  un¬ 
knowns  in  this  formulation.  One  approach  is  to  interpolate  both  unknown 
functions  Y  and  p  in  the  piecewise  bilinear  basis  and  solve  for  vectors  Y  and 
P  from  a  variational  problem  that  incorporates  the  density  relation  f(p)  =  0 
as  a  constraint  condition  by  the  use  of  Lagrange  multipliers.  This  approach 
has  not  been  developed  since  the  earlier  variational  analysis  indicated  that 
the  density  variation  expressions  are  zero. 

The  procedure  advanced  and  applied  here  involves  a  coupled  iteration  for 
position  Y(x,4')»  and  density  pCx,1}1).  The  function  Y (x , y)  is  approximated  on 
the  finite  element  discretization,  and  substituted  in  the  variational  problem 
8  Y 

determining  - r  =  0.  This  leads  to  a  nonlinear  algebraic  system  of  finite 

3Y! 

J 

element  equations  for  Y  and  dependent  on  the  unknown  density  p(x,¥).  If  the 
correct  density  were  known,  then  solution  of  the  algebraic  system  could  be 
achieved  by  a  gradient  method  such  as  Newton-Raphson  iteration.  Starting  with 
an  initial  streamfield  guess,  the  method  iteratively  improves  the  position 
vector,  solving  a  linear  Jacobian  system  at  each  iterative  step.  The  approach 
can  be  extended  to  include  the  unknown  density  as  follows.  Select  an  appropri¬ 
ate  streamline  guess  for  the  solution,  based  on  the  wall  geometry.  The  sca’ar 
density  equation  f (p)  =  0  is  solved  iteratively  with  this  proscribed  stream¬ 
line  position  to  determine  the  corresponding  density  fie’d.  With  this  density 
field  (as  if  it  were  the  exact  result),  the  Newton-Raphson  iteration  for  pos  - 
tion  is  commenced.  Interrupting  the  iteration  at  an  appropriate  point  returns 
an  improved  streamline  solution.  The  density  field  is  recomputed  for  the  new 
streamfield,  and  the  Newton-Raphson  iteration  for  position  recommences  with  the 
new  density  field.  The  coupled  process  of  scalar  density  iteration  and  vector 
position  iteration  is  repeated  to  convergence  within  a  prescribed  error  toler¬ 
ance  on  the  nonlinear  system.  The  linking  of  iterations  is  depicted  in  the 
chart  of  Fig.  3. 

2.  DENSITY  ITERATION 


Given  a  streamline  configuration  Y(x,’f),  the  density  function  satisfies 
the  nonlinear  scalar  equation,  Eq  (29).  This  equation  can  be  solved  at  each 
node  point  (i,j)  where  H,  £  and  Y  are  known  values. 

Newton  iteration  for  the  root  a  of  the  scalar  equation  f  (p)  =  0  is 

fV 

P,  =  P,  -  -  ,  provided  f!(a)  4  0  (48) 

k  f1^) 
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.7ipure  3.  Coupling  oi  Density  and  Streamline 
Calculations  as  an  Iteration 
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Differentiating  Eq  (29) 


f'(P> 


2  (2H  -  ^-)p  -  (Y  +  1)  2h  p 


Y 


y2'H  '  '  '  *■ 7  -‘-Oh 

Substituting  in  Eq  (48)  and  simplifying, 


k+1 


2H 


C2/Y2  2h  0  Y+1  1  +  ** 

Y  "  2h°  Pk  “  (YYy) 2 


2 (2H  -  C  /Yz)p  2-  2 (Y  +  l)h  p  Y+1  J 

k  ok 


(49) 


The  iteration  is  considered  to  have  converged  when  the  difference  in 
successive  iterates  is  less  than  a  specified  error  tolerance.  The  starting 
guess  for  the  initial  calculation  is  P  =  1.  Thereafter,  the  iterates  begin 
from  the  previous  flow  field  estimate.  The  scalar  iteration  is  carried  out  for 
each  node  point  in  the  mesh. 


3.  POSITION  VECTOR  ITERATION 


Let  J  be  the  Jacobian  matrix  J ^  .  The  vector  Newton  iteration 

follows  directly,  as  in  the  scalar  case,  by  Taylor  series  expansion  of 
g (Y  +  AY)  about  Y.  Then  the  vector  iteration  is 


Yk+1 


(50) 


In  practice,  the  Jacobian  derivatives 
difference  approximations, 


3gi/3Y_. 


are  determined  as  finite 


8i(Y»’-"’Yi-i’Yi  +  AV  Vi’  —  ’V  ~  gi(Yt 


. u 


9Y  j 


AY. 
J 


k+1 

The  iterate  Y  is  obtained  as  the  solution  of  the  linear  system 


(51) 


JAY 


k+1 


g(Yk) 


(52) 


J  .  Avk+1  vk+l  ,  3Y 

with  AY  =  Y  -  Y  and  g  =  -rrr 

-  ~  p  Ji 

p  P 

The  coefficient  matrix  J  is  positive  definite,  symmetric,  and  sparse. 
Positive  definiteness  results  from  the  convexity  of  the  functional  v  whose 

stationary  value  is  sought.  Writing  gk  =  E={E™}  for  node  (m,n)  the  symmetry 

follows  from 
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The  sparseness  is  evident  as  the  algebraic  equation  corresponding  to  any 
given  node  has  possible  non-zero  coefficients  only  for  those  nodes  belonging 
to  elements  adjacent  to  the  node  concerned.  Finally,  it  is  important  to 
note  that  a  judicious  numbering  of  the  nodes  will  produce  an  algebrt  problem 
that  is  neatly  banded.  Then  highly  efficient  bandsolvers  can  be  employed  to 

solve  Eq  (52)  for  AY^*1  .  In  the  particular  instance  of  the  throughflow  problem, 
the  characteristic  span  between  hub  and  tip  walls  is  small  compared  with  the 
axial  length  between  the  remote  upstream  and  downstream  stations.  Furthermore, 
the  solution  is  well-behaved  in  the  radial  direction  so  that  a  coarse  mesh 
gradation  in  this  direction  is  appropriate;  that  is,  few  T  lines  are  necessary. 

On  the  other  hand,  there  will  be  a  much  larger  number  of  axial  mesh  lines 
corresponding  to  a  graduated  mesh,  especially  fine  near  the  actuator  discs  and 
becoming  coarse  towards  the  remote  boundaries.  With  this  knowledge  of  the  flow 
behavior  and  mesh  natuie,  the  optimal  nodal-numbering  scheme  is  evident: 
beginning  at  the  remote  upstream  boundary,  number  radially  across  the  stream¬ 
lines.  Then  node  k,  cr  rresponding  to  mesh  pair  (m,n)  for  a  mesh  of  N  T-lines 
and  M  axial  mesh  lines  (with  M  »  N),  is  (m-1)  N  +  n.  The  half-bandwidth  is  N 
and  the  coefficient  matrix  is  block-tridiagonal .  Only  the  bandedness  will  be 
exploited  in  this  algorithm  as  there  is  non-zero  "fill"  between  the  tridiago¬ 
nal  blocks  during  decomposition. 


For  direct  solution  of  symmetric,  positive-definite  systems,  the  Cholesky 
method  is  preferable.  The  Cholesky  decor.oosition  is  the  triangular  decomposi¬ 
tion,  J  =  L  Lt,  where  L  is  lower  triangui  r.  In  the  program  code  an  equivalent 
form  of  Cholesky  decomposition  was  used,  expressing  J  as  (JTy  where  IJ  is  upper 
triangular.  Thio  allowed  more  convenient  computation  and  display  of  the  com¬ 
pactly  stored  jacobian  and  its  decomposed  form.  Writing  J  =  U^U  so  that 

UTU  AYk+1  =  g(Yk),  implies: 

UTS  = 

W+l 

and  UAY  =  6 


forward  reduction 
back  substitution, 


to  determine  the  solution 


Only  the  non-zero  upper  half-band  of  J  need  be  stored  and  U  is  overstored 
on  J  as  decomposition  proceeds. 


4.  STABILITY  AND  CONVERGENCE 


Fixed  point  theory  for  scalar  and  vector  iterations  is  quite  extensively 
developed  and  the  Newton-Raphson  iteration  has  been  investigated.  Standard 
contraction  arguments (*)  indicate  that  the  iteration  converges  if  the  initial 
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guess  is  within  the  contraction  interval.  The  streamline  position  iteration  is 
generally  well-behaved,  as  the  wall  geometry  suggests  a  reasonable  starting 
guess  which  is  coded  directly  in  the  program  and  need  not  be  developed  as  data. 
The  density  Iteration  requires  more  care.  Given  an  initial  streamline  configura¬ 
tion,  the  starting  guess  is  incompressibility,  p  =  1,  and  Newton  iteration  of 
the  nonlinear  scalar  equation  at  each  grid  point  is  generally  convergent. 

However,  in  regions  where  there  is  a  dramatic  change  in  density,  such  as  on  thf 
hub  axis  immediately  beyond  a  stator  imparting  large  swirl,  the  density  itera¬ 
tion  may  not  converge.  An  inadequate  streamline  configuration  estimate, 
together  with  a  poor  initial  density  iterate,  may  predict  meridional  Mach  num¬ 
ber  values  greater  than  unity  when  the  real  physical  flow  has  a  subsonic  meri¬ 
dional  Mach  number. 

To  stabilize  the  density  iteration,  the  successive  Newton  iterates  are 
scrutinized:  if  the  iteration  at  a  particular  grid  point  begins  to  diverge, 

the  physical  constraint  to  subsonic  meridional  Mach  numbers  is  enforced  by 
setting  axial  velocity  W  =  1  and  radial  velocity  U  =  0.  This  determines  a 
reasonable  interim  approximation  at  the  troublesome  point.  Subsequent  posi¬ 
tion  iterations  improve  the  streamline  guess,  and  the  later  density  iteration 
generally  converges,  with  the  exceptions  of  borderline  cases  where  the  physical 
flow  is  very  demanding  and  sonic  meridional  Mach  numbers  are  being  closely  ap¬ 
proached.  T'iis  approximate  measure  then  serves  to  stabilize  the  numerical 
iterations  and  extends  the  numerical  stability  to  include  essentially  the  full 
range  of  flows  of  interest. 

The  graphs  in  Figs.  4  and  5  plot  the  error  at  a  point  on  the  hub  wall 
beyond  a  stator,  for  several  computations  at  various  incident,  free-stream  Mach 
numbers.  The  graphs  show  very  clearly  the  influence  of  the  stabilizing  density 
approximation. 

Another  technique  for  stabilizing  the  density  calculation  is  to  moderate 
the  flow  transition  through  the  blade  row.  This  can  be  done  readily  by  using 
several  actuator  discs  closely  nested  with  a  very  fine  separating  grid,  rather 
than  use  a  single  disc.  Such  an  example  is  later  examined  in  SECTION  VI. 

The  exploration  of  the  density  calculation  brings  forth  the  question  of 
stability  of  the  global  coupled  iteration  process,  rather  than  that  of  the  posi¬ 
tion  or  density  iterations  alone.  In  particular,  should  the  position  iteration 
at  a  fixed  density  be  iterated  to  convergence  (a  quasi-incompressible  analysis) 
and  the  density  then  updated,  or  should  the  position  iteration  be  terminated 
earlier  Early  computations  explored  the  former  scheme  and  it  was  found  that 
in  many  cases  the  coupled  process  diverged.  The  use  of  under-relaxation  at¬ 
tempting  complete  position  iteration  was  similarly  unsuccessful  in  these  cases. 
Limiting  the  number  of  position  iterations  and  successively  reducing  the  limit 
eventually  led  to  a  coupled  sequence  that  was  stable  and  converged.  The  best 
results  were  obtained  using  a  single  position  iteration  to  improve  the  stream¬ 
line  pattern.  This  analysis  indicated  the  necessary  close  coupling  of  the 
numerical  schemes. 

Computationally,  the  global  rate  of  convergence  of  the  coupled  calcula¬ 
tions  is  only  linear,  whereas  the  Newton  processes  are  asymptotically  quadra¬ 
tic.  For  each  case  the  rate  of  convergence  of  both  position  and  density  is 
somewhat  slower  for  the  first  few  iterations.  The  density  point  computations 
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NUMBER  OF  COMPLETE  COUPLED  CYCLES 


NUMBER  OF  COMPLETE  COUPLED  CYCLES 

Figure  5.  Plots  of  Peak  Meridional  Mach  Number  on  the  Hub  Streamline  Immediately  Following  the  Stator 


required  several  iterations  in  the  starting  calculations.  However,  for  most 
incident  Mach  numbers  and  away  from  heavily  loaded  discs,  subsequent  point 
density  calculations  required  only  one  or  two  iterations  to  give  six  signifi¬ 
cant  figure  accuracy. 

Finally,  the  scheme  is  naturally  sensitive  to  the  incident  Mach  number. 

As  Mach  number  increases  the  calculation  becomes  more  demanding,  stabilization 
is  necessary,  and  the  rate  of  convergence  decreases.  This  is  due  in  part  to 
the  physical  nature  of  the  higher  Mach  number  flows,  attainment  of  a  sonic 
meridional  Mach  number  being  associated  with  "choking."  In  addition,  sequences 
of  flows  computed  at  successively  higher  incident  Mach  numbers,  on  the  same 
configuration  and  at  the  same  loading,  all  used  the  same  starting  streamline 
pattern  This  pattern  is  best  for  nearly  incompressible  flows  and  deteriorates 
with  increased  incident  Mach  number.  Thus  it  takes  increasingly  more  iterations 
to  improve  the  starting  guess.  Relaxation  and  extrapolation  techniques  are 
possible  ways  to  accelerate  convergence  but  have  not  yet  been  explored,  the 
present  scheme  affording  adequate  results  in  fewer  than  ten  coupled  iterations 
for  most  applications. 
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SECTION  V 

THE  THROUGHFLOW  PROGRAM 


1.  BASIC  PROGRAM  STRUCTURE 

The  throughflow  problem  is  specified  by  input  data  to  the  program,  the 
finite  element  idealization  is  established,  and,  beginning  with  a  preset  stream- 
field,  the  coupled  iterations  in  density  and  position  are  carried  out.  Each 
position  iteration  is  associated  with  solution  of  a  banded  linear  system  of 
equations.  The  general  computational  procedure  and  subroutine  inter-relations 
are  depicted  in  the  macrochart  of  Fig.  6. 

The  function  of  each  module  and  its  relation  to  other  stages  in  the  layout 
are  described  briefly  below. 

a.  Input 

The  input  instructions  reside  within  the  controlling  program,  TCOMFLO. 
Input  data  defining  the  problem  is  read  from  cards  for  each  flow  problem 
analyzed. 

b .  TCOMFLO 

This  is  the  main  calling  program.  It  determines  the  finite  element 
idealization  in  the  (x.'F)  coordinate  system,  and  controls  the  data  input  and 
solution  output  steps. 

C.  GUESS 

The  initial  guess  for  YCx.'f)  at  the  mesh  points  is  prescribed. 

d.  RSOLVE 

This  contains  the  density  iteration  and  position  vector  iteration, 
returning  the  final  values  to  the  main  controlling  module,  TCOMFLO. 

e.  BANSOL 

This  subprogram  is  called  from  RSOLVE  and  bandsolves  the  linear 
Jacobian  system  for  AY. 

f .  MERMACH 

The  meridional  Mach  number  is  computed  at  each  node  point  between 
streamfield  position  iterations  to  monitor  convergence  of  the  coupled  itera¬ 
tions. 


g.  Output 

The  final  streamline  coordinates  and  associated  nodal  velocity  com¬ 
ponents  are  printed  by  axial  mesh  section. 

Important  details  of  the  program  are  now  developed. 
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Input 


Output 


Figure  6. 


Modular  Subroutine  Layout  of  Programs 
for  Throughflow  Calculations 


2. 


PROGRAM  METHOD 


F’ 


1 


Finite  element  mesh  generation.  The  sections  of  uniform  flow  between  discs 
are  delimited  by  variables  set  in  the  program  code.  For  a  single  stage,  the 
separation  between  stator  and  rotor  is  input,  and  a  uniform  fine  mesh  established 
between  the  disc  pair.  Away  from  the  discs,  the  grid  spacing  increases  with 

■  K  Xj  where  AXj+^  “  xj+i  “  xj  downstream  and  Ax^_^  =  xj  ~  xj  1  uPstream 

and  <  is  set  in  the  program  as  1.25.  Most  of  the  results  described  later  em¬ 
ployed  this  restricted  mesh  algorithm.  For  general  purposes,  problems  with 
several  stages  will  require  more  flexible  data-generation  schemes  or  simply  in¬ 
put  the  values  (x  }  as  in  Ref.  3. 


Initial  Streamfield  Pattern.  The  starting  pattern  is  determined  from  the 
hub  and  tip  wall  shapes  and  the  uniform  flow  at  the  remote  upstream  boundary. 
At  the  spstream  boundary  the  flow  is  uniform  with  axial  velocity  W_  =  1  and 

1  3*^  m 

f  =  1.  Then  the  axial  velocity  satisfies  pW  =  -  —  -r—  and  leads  in  the  (x  =  'F) 

1  2  Y  1  ? 

system,  to  — (Y  )  =  1.  Integrating  with  respect  to  t  gives  —  Y  =  -f  so  that 


Y(x0,  V) 


(53) 


and  determines  the  f-line  distribution  at  the  remote  upstream  boundary.  The 
'{'-line  spacing  for  the  initial  guess  is  chosen  so  that  the  separation  of  adja¬ 
cent  streamlines  is  constant  along  any  radial  mesh  line  from  hub  to  tip.  In 

the  (x,T)  plane,  then  -  4^  =  -  -fj^i+i]2  “  (Yi)2]'  Clearly  the  spacing 

AY  will  vary  from  one  mesh  line  j  to  the  next  as  the  hub  wall  f(x)  or  tip  wall 
g(x)  change.  Then  for  a  mesh  line  at  x^  , 

,  r—  (g(x  )  -  f(x  )) 

Y(xj>V  -  Y{  =  f(xj)  +  (V^  -  ^xo»  (g(xJo)  --fO^T  (5A) 


The  Density  Iteration. 


The  density  at  any  node  point  (x^,^)  is  deter¬ 


mined  by  Newton  iteration  applied  to  the  equation  f(P)  =  0  as  described  in  the 
preceding  section.  For  highly  rotating  flows  at  high  incident  Mach  number, 
the  density  iteration  may  require  stabilization.  This  is  particularly  true  in 
the  early  stages  of  coupled  iteration  as  the  streamline  pattern  will  be  poor 
in  some  regions.  The  density  calculation  is  indicated  in  the  chart  of  Fig.  7. 


Streamline  Position  Adjustment.  The  variational  finite  element  formula¬ 
tion  for  strealine  position  Y(x,40  leads  to  the  nonlinear  system  of  equations 


3V_ 

K 

8  y  _ 2. 

5yj 

m,n  9Y1 

j 

0 
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Writing  this  as  g(Y)  -  0,  the  solution  Y  is  sought  by  Newton-Raphson  iteration. 

V.  *  R  k. 

Let  Y  correspond  to  the  k  th  intermediate  solution  pattern.  Then  g(Y  )  =  e 

+  Q,  and  is  the  error  vector.  The  solution  adjustment  is  achieved  as  a  single 

step  of  the  Jacobian  system  iteration 


JAY 


k+1 


-e 


In  SECTION  IV,  the  element  contributions  to  3V/3Y^  were  developed  and  per¬ 
mutation  relations  noted.  These  permutation  relations  allow  an  elegant  coding 

k  k 

for  accumulation  of  g  =  g  at  each  node  point,  and  similarly  for  the  Jacobian 

evaluation.  A  single  statement  function  definition,  FCY™  ,  Y1",,  ,  ,  Y111"^) 

n  n+1  n  n+1 

determines  the  element  contribution  to  a  given  node,  and  permutation  of  the 

Ct 

corner  coordinates  (YR)  ascertain  the  other  nodal  contributions  from  the  element 
(Fig.  8).  P 


The  compuational  procedure  is  outlined  in  Fig.  9.  Using  the  density  field 
calculated  from  the  current  streamline  position,  the  new  adjusted  streamline 
pattern  is  computed.  Proceeding  sequentially  by  elements,  the  element  contri¬ 
bution  of  8V™/3Yj  to  each  of  the  four  corner  nodes  is  determined  by  use  of 

the  statement  function.  The  nodal  contributions  are  accumulated  to  give  the 
error  vector 


/  3V  1 

1 3Yj/ 


-  g or>  - 


and  the  infinity  norm  of  e  is  tested  to  determine  if  the  solution  has  con¬ 


verged.  The  Jacobian  J 


3g  /3Y  is  formed  numerically  by  differencing, 
pq  P  q 


3^  g(Y„  ...Yp.!  ,  Yp  +  AYp>  Yp+1 . Y^  -  g^,..^) 

3Y  ~  AY 

q  p 


Again,  this  difference  quotient  is  formed  sequentially  by  elements,  evaluating 

g (Y  , . . . , Y  +  Ay  ,...,Y  )  and  g(Y  , . . . ,Y  )  by  accumulating  element  contributions 
lppN  1  N  k  _g 

to  the  nodes  as  in  the  preceding  calculation  of  e  .  The  value  AY^  =  10  is  set 

in  the  program  code.  Only  the  banded  upper-half  of  the  jacobian  is  stored 
during  formation.  The  linear  jacobian  system  is  solved  by  banded  Cholesky 
decomposition,  followed  by  forward  and  backward  substitution  sweeps  on  the 

k 

right-hand  side  vector  g  .  During  decomposition,  only  the  upper  triangular 
decomposed  matrix  is  needed,  and  is  stored  over  the  coefficient  matrix  J  as 

k+1 

decomposition  proceeds.  Solution  for  AY  determines  the  adjusted  streamline 
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pattern.  After  the  new  position  solution,  we  return  to  the  density  calculation 
and  the  next  coupled  iteration  commences. 


3.  PROGRAM  USAGE 


Idealization  Description.  The  turbine  configuration  is  described  in  two 
dimensions,  and  consists  of  a  sequence  of  sections  of  continuous  flow  separated 
by  line  interfaces  representing  actuator  discs.  The  lower  and  upper  boundaries 
correspond  to  hub  and  tip  surfaces  respectively,  and  the  remote  left  and  right 
boundaries  simulate  uniform  flow  regions  far  upstream  and  far  downstream.  The 
flow  is  characterized  by  a  discrete  streamline  pattern,  and  the  problem 
formulated  in  the  (x , 4')  plane  to  solve  for  streamline  location.  In  (x,T) 
coordinates  the  domain  is  rectangular  and  a  finite  element  idealization  of  rec¬ 
tangular  is  introduced,  with  grid  lines  4*  =  4*  determined  in  the  program,  and 

x  =  x.  -*nt)ut  as  data.  Angular  momentum  and  enthalpy  far  upstream  and  at  the  discs 

are  prescribed  for  each  section  of  continuous  flow.  A  typical  configuration  with 
three  discs  and  varyinp  hub  and  tip  walls  is  shown  in  Fig.  10.  The  finite  ele¬ 
ment  idealization  in  the  (x,T)  plane  is  indicated  in  Fig.  11. 


A  typical  rectangular  element,  shaded  in  section  2  of  Fig.  11,  is  defined 
by  enclosing  horizontal  Y-lines  and  vertical  x-lines.  The  section  limits  Nz^  ^ 
with  £  =  1,...  sections  and  k  =  1,2  delimit  each  domain  of  continuous  flow  ’ 
and  are  the  x-l’.ne  index  pairs  marked. 


Computer  application  requires  specification  of:  the  idealization  as  g’-id 
point  ordinates  in  the  (x,Y)  plane  and  identification  of  flow  sections;  angular 
momentum  and  enthalpy  at  each  interface  between  adjacent  sections.  The  rectangu¬ 
lar  domain  in  the  (x,'f)  plane  simplifies  the  finite  element  representation  and 
for  flexibility  the  prescription  of  the  x-partition  is  required  as  card  input. 

The  ’f-partition  is  set  in  the  program  and  corresponds  to  AY  =  constant  at  the 
upstream  boundary.  The  fully-inf ini te  pioblem  is  modelled  as  a  finite  but 
extensive  region,  and  natural  boundary  coi  ditions  at  the  remote  ends  of  the  flow 
approximate  the  free-stream  conditions. 


List  of  Program  Variables. 


ALFA(K):  Coefficient  of  momentum  fit  across  disc,  yV  =  -OtT  +  S 


BETA(K) : 
C(N,M): 

E (N) : 
F(J): 
G(J): 

H (I ,  K) : 

I  DEN: 
ITMX: 
MNMAX: 
NR: 


Coefficient  of  momentum  fit  across  disc,  vV  =  -«4'  + 

Upper  half-band  of  symmetric  jacobian  matrix 

Error  residual  vector  of  nonlinear  system 

Vector  defining  lower  wall  (hub)  ordinate  boundary,  f (x) 

Vector  defining  upper  wall  (tip)  ordinate  boundary,  g(x) 

Total  enthalpy  on  streamline  I  in  section  K 

Maximum  number  of  cycles  of  coupled  calculation 

Maximum  number  of  position  iterations  per  cycle,  usually  1 

Maximum  number  of  density  iterations  per  cycle 

Row  dimension  of  jacobian,  C 


t 


Indices : 


K  =  section,  I 


streamline,  J  =  x-line,  N  =  node  (I,J) 
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Tip 


Disc 


Figure  10.  Multiple  Disc  Throughflow  Configuration 
Typical  Rectangular  Element 


Figure  11.  Idealization  of  Turbine  Flow  Problem  to  a  Finite  Element 
Representation  in  the  (x,40  Plane 


NS: 

NY: 

NZ(K,L) : 
NZMX: 
OMEGA (K) : 
RON ( I , J ) : 
T(I,K): 

X(I, J) : 

Y  (I) : 
Z(J): 


Number  of  flow  sections 
Number  of  horizontal  'f-lines 
Section  delimiters  in  pairs,  L  =  1,2 
Number  of  vertical  x-lines 

Angular  velocity  introduced  to  each  section  (T^) 
Density  at  node  (I,J) 

Section  angular  momentum  on  each  T-line,  T  =  (yV)‘ 

Solution  array  of  f-line  ordinates,  (y|) 

Stream  function  values,  ¥ 
x-grid  abscissas,  x 


c2 


Input  Preparation.  The  input  scheme  is  described  sequentially  by  card 
type.  Numbers  in  parentheses  following  each  listed  number  indicate  the  number 
of  card9  of  this  type  required.  Variables  are  listed  within  each  card  type  in 
order,  and  correspond  to  the  column  fields  as  sequenced. 


(a) 


Card  Type  1 


(1) 


Cols : 

Format: 

Variables: 

Description: 


Comments: 


1-5,  6-10,  11-15,  16-20 
415 

NY,  NZMX,  NS,  NF 

NY  -  Number  of  T-lines 
NZMX  -  Number  of  x-lines 
NS  -  Number  of  sections 
NF  -  Flag:  0  ■  enclosed,  1  =  free  exit 
Use  NF  ■  0  for  this  program  version. 


(b) 


Card  Type  2 

Cols: 

Format: 

Variables: 

Description: 


Comments: 


(NZMX) 

1-5,  6-15,  16-25,  26-35 
15,  3F10.5 
J,  Z(J),  F (J) ,  G (J) 

J  -  Vertical  grid-line  z-index 
Z(J)  -  Z-abscissa 
F(J)  -  Hub-surface  boundary 
G(J)  -  Tip-surface  boundary 

Repeat  each  card  to  end  of  z-lines;  J  *  1,  NZMX 


(c)  Card  Type  3 


([hS/3]  +  1) 


Cols: 

Format: 

Variables: 

Description: 

Comments: 


1-10,  11-20,  21-30,  31-40,  41-50,  51-60 
3(2F10.5) 

NZ(L,K) 

NZ(L,K)  -  Flow  section  delimiters  in  pairs, 

L  -  1,2 

Set  up  in  sequence  with  3  sections  per  card  till 
end  of  sections,  e.g.,  the  first  card  will  have 
in  order:  NZ(1,1),  NZ(2,1);  NZ(1,2),  NZ(2,2); 
NZ(1, 3) ,  NZ (2 , 3)  . 
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(d) 


(e) 


(f) 


(g) 


(h) 


Card  Type  4 

(1) 

Cols: 

1-2 

Format : 

12 

Variables : 

I  FLAG 

Description: 

(=  1  implies  read  card  types  6b  and  7  ) 

L  a®  (j4  1  implies  read  card  type  6a  or  ly  ) 

Comments : 

IFlag  determines  mode  of  prescribed  enthalpy  and 
momentum  within  sections. 

Card  Type  5 

(1) 

Cols: 

1-10 

Format : 

F10.5 

Variables : 

AMD 

Description: 

AMD  -  Incident  Mach  number  for  upstream. 

Card  Type  6a 

(NS) 

Cols : 

1-10,  11-20,  21-30 

Format : 

3F10.5 

Variables : 

ALFA(K) ,  BETA(K) ,  OMEGA (K) 

Description: 

Momentum  (rv)  =  -ALFA(K)*T  +  BETA(K) 

ALFA(K)  —  Swirl  component 

BETA(K)  -  Translation  component 

OMEGA(K)  -  Angular  velocity 

Comments : 

Linear  fit  to  momentum  introduced  by  disc  at  sec 
tion  entrance.  Repeat  each  card  to  end  of  sec¬ 
tions  K  =  1,  NS.  Repeat  sequential  sets  for  para 
metric  studies  (see  Recommendations). 

Card  Type  6b 

([NS/8]  +  1) 

Cols : 

1-10,  11-20,  21-30,  31-40,  41-50,  51-60,  61-70, 
71-80 

Format : 

8F10.5 

Variables : 

OMEGA (K) 

Description: 

OMEGA(K)  -  Angular  velocity  introduced  to  each 
flow  section 

Comments : 

Repeat  to  end  of  sections. 

Card  Type  7 

( [NS*NY /8] ) 

Cols : 

1-10,  11-20,  21-30,  31-40,  41-50,  51-60,  61-70, 
71-80 

Format : 

8F10.5 

Variables : 

T(I,K) 

Description : 

T(I,K)  -  Angular  momentum  (rv)  at  T-line  I  for 
section  K 

Comments : 

Repeat  each  card  to  end  of  sections,  K  =  1,  NS 
in  sequential  order  by  streamlines  1=1,  NY. 
Repeat  sets  6b  and  7  for  parametric  studies  (see 
Recommendations)  . 
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Recommendat  Ions .  Multiple  data  sets  can  be  processed  and  are  to  be 
separated  by  END- OF-RECORD  cards  (7-8-9  multipunch  in  column  1).  Parametric 
studies  for  varied  section  angular  velocities  and  momenta  can  be  run  within 
each  full  data  set  by  duplicating  card  types  6  and  7.  Also,  the  incompressible 
flow  calculation  is  available  as  a  default,  by  setting  MNMAX  =  0  as  the  itera¬ 
tion  limit  for  density  calculations. 

Restrictions .  Arrays  currently  are  dimensioned  in  the  main  program 
TCOMFLO  to  accommodate  at  most:  15  streamlines,  42  z-lines,  3  sections.  Note 
that  NR  is  set  to  630,  the  current  row  dimension  of  C(M,N).  The  corresponding 
dimension  statements  can  be  altered  readily.  Subprograms  are  variably  dimen¬ 
sioned  and  require  no  adjustment.  Array  usage  is: 

NZ  (2 ,  i  NS),  F(i  NZMX) ,  G(i  NZMX)  ,  H(NY ,  *  NS),  T  (NY ,  i  NS), 

X(NY ,  £  NZMX),  Y(>  NY),  Z(£  NZMX),  C(NR,  £  NY  +  2),  E(>  NODE), 

ALFA(£  NS),  BETA(^  NS),  OMEGA (^  NS),  RON (NY ,  *  NZMX) 


To  change  the  array  dimensions  the  relevant  source  cards  are:  cards  2, 

3,  4  and  12  in  program  TCOMFLO.  The  iterations  are  controled  by  the  following 
variables  set  in  RSOLVE: 


IDEN  =  10 

> 

the  maximum  number  of  coupled  cycles 
allowed 

MNMAX  =  5 

» 

the  maximum  number  of  Newton  itera¬ 
tions  for  density  at  each  grid  point 

ITMX  =  1 

y 

the  maximum  number  of  position  itera¬ 
tions  of  the  finite  element  system 

ETOL  =  10~5 

_  c 

y 

error-bound  tolerance  for  iterate 
convergence 

DX  =10 

y 

difference  increment  for  jacobian  de¬ 
rivative 

Convergence  is  accepted  if  the  infinity  norm  of  E  is  less  than  ETOL.  If 
more  than  IDEN  cycles  are  necessary  the  program  terminates  with  a  message,  and 
prints:  the  current  contents  of  the  streamline  position  solution,  the  last 

position  adjustment,  and  the  velocity  components  at  each  node. 

Diagnostics .  The  following  diagnostic  comments  may  appear  during  program 
execution: 


"IMPROPER  ARRAY  SIZE  FOR  EXECUTION"  —  accompanied  by  the  input 
NR,  NS,  NY,  NZMX,  NODE  values.  This  is  triggered  by  either  in¬ 
correct  data  (prescribing  NR)  or  the  need  to  redimension  arrays 
in  the  program  to  accommodate  a  "larger"  problem.  The  program 
terminates  after  printing  the  diagnostic. 

"SOLUTION  DID  NOT  CONVERGE"  —  the  maximum  number  of  full 
cycles  was  completed  and  ||fi|  had  not  yet  decreased  to  the 
prescribed  tolerance  ETOL.  The  control  returns  to  the  calling 
program,  whereupon  the  velocity  field  is  calculated  and  printed 
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"ENCOUNTERED  MERIDIONAL  MACH  NUMBER  IN  EXCESS  OF  UNITY"  — 
printed  from  subroutine  MERMACH,  this  indicates  that  the  cur¬ 
rent  position  and  density  is  predicting  a  supersonic  meridional 
Mach  number.  For  physically  realistic  problems  this  is  cor¬ 
rected  in  later  cycles  by  density  and  position  adjustment.  Com¬ 
putation  continues  to  complete  the  analysis.  Physically  ill- 
posed  problems  will  lead  to  an  arithmetic  error  in  the  band- 
solver,  caused  by  divergence  of  the  solution. 

"DENSITY  ADJUSTED  TO  STABILIZE  COUPLED  PROBLEM"  —  this  occurs 
when  the  position  guess  is  so  poor  that  locally  the  Newton 
iteration  for  density  fails.  A  physical  argument  resets  the 
density  to  an  approximate  estimate  and  calculation  proceeds. 


Subprograms  and  Calling  Sequences. 


GUESS:  Call  GUESS  (X,F,G,NZMX,NY,Y,XT,XH,YT,YH}  [also  call 

COORD  (X,F,G,NZMX,NY,Y,XT,XH,YT,YH)  as  a  second 
entry  point].  GUESS  determines  the  V-line  distri¬ 
bution  and  the  starting  streamfield  pattern. 

RSOLVE :  Call  RSOLVE  (NF,NR,NS,NY,NZ,C,E,F,G,H,T,PE,X,Y,Z, 

RON).  RSOLVE  solves  for  density  and  position  adjust¬ 
ments  as  an  iterative  cycle. 

MERMACH:  Call  MERMACH  (NS ,NR,NZ,H,T ,NY,X, Y ,RON,Z,NDEN)  . 

MERMACH  calculates  the  meridional  Mach  number  at  each 
node  in  the  idealization. 


FANSOL:  Call  BANSOL  (NBAN,NODE,NR,C,E) .  BANSOL  solves  the 

symmetric  positive-definite  Jacobian  system  for  each 
position  adjustment,  using  a  banded  Cholesky  decom¬ 
position  and  substitution  sweeps. 


Machine  Requirements.  The  program  has  been  developed  on  the  CDC  6400 
computer  at  the  University  of  Washington.  It  is  coded  completely  in  FORTRAN 
with  no  machine-dependent  instructions,  and  so  can  be  adapted  to  alternate 
hardware.  With  current  array  dimei  ions  the  program  requires  42,000  octal 
central  memory  storage  locations  and  takes  approximately  9u  seconds  CDC  6400 
central  processor  time  for  solution  at  the  maximum  array  size  allowed.  A 
detailed  breakdown  of  timing  for  sample  calculations  follows. 

Timing .  Timing  profiles  for  a  sample  problem  with  coarse  and  fine  mesh 
idealizations  are  shown  in  Figs  12  and  13.  The  calculation  is  broken  down  to 
its  main  components  —  the  density  iteration  (D) ,  jacobian  calculation  (J), 
bandsolver  (B) ,  position  adjustment  (P) ,  and  total  coupled  iteration  (T) .  The 
graphs  indicate  that  the  density  iteration  is  quite  inexpensive  and  confirm  the 
earlier  reasoning  in  formulating  the  coupled  process. 

Output .  An  echo-check  of  the  input  data  is  displayed;  the  density,  posi¬ 
tion,  and  error  arrays  are  printed  during  each  cycle  together  with  the  peak 
meridional  Mach  number  on  each  T-line  and  in  each  section.  Finally,  the  posi¬ 
tion,  stream  function,  axial  velocity  and  tangential  velocity  are  printed  by 
f-line  beginning  at  the  hub  and  for  each  x-line  segment. 
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CENTRAL  PROCESSOR  TIME, 


r  r 


R  -  10/7 
Parallel  Walls 
S-R  separation  0.14 
a  “  0 
6  -  6.66 
•  0.182 


T:  Total  (P+D) 


Position  Adjustment 
(Includes  J  and  B) 


J:  Jacobian  Calculation 


B:  Bandsolver 


D:  Density  Iteration 


2  4  6  8  10 

NUMBER  OF  COMPLETE  COUPLED  CYCLES 


Figure  12.  Timing  in  Central  Processor  Seconds  for  Coarse  Mesh  (8  x  24) 
Calculation 


NUMBER  OF  COMPLETE  COUPLED  CYCLES 


Figure  13.  Timing  in  Central  Processor  Seconds  for  Fine  Mesh  (8  x  42) 
Calculation 
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Deck  Structure  for  Compile  and  Execute  from  fource  Code. 


JOB  CARD 

RUN (S) 

(Compile) 

Job  Control 

LGO 

END-OF-RECORD  CARD 

(Load  and  execute) 

Record 

PROGRAM  CODE 
END-OF-RECORD  CARD 

1st  Data  Set 

END  OF  RECORD  CARD 

Source  Deck 

FINAL  DATA  SET 
END-OF-RECORD  CARD 
END-OF-FILE  CARD 


Sample  Problem.  A  simple  throughflow  example  with  single  stage  and  varying 
tip  surface  is  identified  in  three  sections  in  Fig  14.  The  corresponding  finite 
element  idealization  in  the  (x.'F)  plane  is  given  in  Fig  15. 


Input . 


8 

24  3 

0 

1 

-9.42 

1.0 

3.0 

2 

-7.04 

1.0 

3.0 

3 

-5.13 

1.0 

3.0 

4 

-3.60 

1.0 

3.0 

5 

-2.38 

1.0 

3.0 

6 

-1.41 

1.0 

3.0 

7 

-0.63 

1.0 

3.0017 

8 

-0.0 

1.0 

3.0085 

9 

0.5 

1.0 

3.0172 

10 

1.0 

1.0 

3.0283 

11 

1.5 

1.0 

3.0410 

12 

2.0 

1.0 

3.0546 

13 

2.5 

1.0 

3.0681 

14 

3.0 

1.0 

3.0808 

15 

3.5 

1.0 

3.0919 

16 

4.0 

1.0 

3.1006 

17 

4.63 

1.0 

3.1074 

18 

5.41 

1.0 

3.1086 

19 

6.38 

1.0 

3.1086 

20 

7.60 

1.0 

3.1086 

21 

9.13 

1.0 

3.1086 

22 

11.04 

1.0 

3.1086 

23 

13.42 

1.0 

3.1086 

24 

16.40 

1.0 

3.1086 

1 

8 

9 

22 

.2 

0. 

0. 

.0 

2. 

0. 

.0 

0. 

0. 
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SUBROUTINE  RSOL VE ( NF ,NR, NS, NY ,NZ , C,E , F, G, H ,T,PE , X, Y , Z , RONI 

SOLVES  A  NONLINEAR  SET  OF  FINITE  ELEMENT  EQUATIONS  GIVING 
A  THROUGHFLOH  MODEL  CF  AN  ANNULAR  TURBCMACHINE  INCLUDING 
VARIABLE  WALL  GEOMETRY  AND  THE  POSSIBILITY  OF  A  NEARBY  EXIT 
WITH  FREE-STPEAMLINES.  THE  TREATMENT  IS  JASED  ON  PAGTUS(X) 
AS  A  FUNCTION  OF  PSI(Y)  AND  AXI AL-DOSIT ION f Z> .  THE  MACHINE 
IS  DIVIDED  INTO  NS  SEGMENTS  OF  CONTINUOUS  FLOW  KITH  NS-1 
ACTUATOR  DISCS  BET  KEEN  INTERNAL  DIVISIONS.  NY  PSI-I.INES 
ARE  IN  EACH  SEGMENT,  AND  NZ(JE.JS)  DEFINES  THE  END  Z-OIVI- 
S  ION  NUMBERS.  WHERE  SPECIFIED,  THE  INNER  WALL  SHAPE  IS 
GIVEN  BY  F ( J )  AND  THE  OUTER  WALL  3Y  G(J>. 

NF=0  EXIT  FAR  COKNSTRtAM  (NO  FREE -STREAMLINES I 
NF=1  LAST  SEGMENT  IS  FREE  (F  A NO  G  TO  BE  DETERMINED) 

H  ( I, JS)  AND  T ( I , JS )  CEFINE  THE  STAGNATION  ENTHALPY  ANO 
TANGENTIAL  VELOCITY  FUNCTIONS  IN  THE  INTERVAL  Y(1),Y(I»-1> 
FOP.  EACH  SEGMENT  ( JS=1,  . .  .  , NS)  .  PE  DEFINES  THE  EXTERNAL 
CONDITIONS  ON  THE  FREE- STREAMLINES.  IT  NEED  NOT  BE  SPECI¬ 
FIED  IF  NF=g  ...  PE=  t  IS  THEN  AUTOMATICALLY  SET. 

NR  IS  THE  ROW  DIMENSION  CF  THE  MATRIX  C(K.L)  AS  SFT  IN 
THF  MAIN  PROGRAM.  C  IS  USEO  TOR  COMPACT  STORAGE  OF  THE 
3ANOCO.  SYMMETRIC  JACOBIAN  MATRIX  USED  IN  SOLVING  THE  NON¬ 
LINEAR  SYSTEM  BY  NEWTON-P.APHSON  ITERATION.  THE  ERROR  RESI¬ 
DUALS  ARE  STORED  IN  E ( K)  BY  PSI-LINE. 

VERSION  ....  FINITE  DIFFERENCE  JACOBIAN  EVALUATION 

DIMENSION  NZ(2,D,F(1),G(1),H{NY;1),T<NY,1) 

DIMENSION  X( 1) ,Y(1) ,Z(1) ,C(NR,1) ,E(1) 

DIMENSION  RON ( NY , 1 ) 

DIMENSION  V ( 4) 

LOGICAL  NTST 
REAL  MMN 

ELEMENT -CHARACTERISTIC  FUNCTION  (MEAN-VALUE  APPROXIMATION) 
0(R1,R2,R3,R4)  =  OZ*(  2*  (  FF-PtJ)  *R4-C  .5*TT*<  1/R3  +  R3/  <  R4*R4)  I 

*  +DY*DY* (  (R4-R2) *ll/(R4* (R4-R3I )  +1/(k2*  C;2-R1) ) )/ (CZ'DZ) 

*  -3.5*(  (1«((R4-R2)/DZ)**2)*(2*R4-K3)/(R4*RA) 

*  * ( 1 ♦ ( (R3-R1) /OZ) **2>/R3>/ (R4-R3) **2> /RN2) *RN 

NY  1-NY-  1  $  NZMX-NZ(2,NSI  »-l 

NBA N=NY »2  T,  NODE=Nv*MZMX 

ITMX=1  3  ETOL=l. t-G  S  0X=i.E-5  t  IDEN=10 

GAM  =  1  ,4 

RF=1.0 

NPRT  =  5 

MNH AX=5 

ITMX=1 

NPSS-2 

I0BN=1C 

IP=ITmx-i 
IF(NF.EO.P)  P£=C. 

IF(NR.GE.NOOE)  GO  TO  1 
WRITE (6,201)  ND»NS,NY»NZMX,NOOE 
CALL  EXIT 
1  CONTINUE 


RUN 
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163 

170 

171 

204 

205 
221 
226 
232 

234 

235 

236 

240 

241 
246 
252 
255 
262 
264 
270 
273 
275 
302 
314 
324 
333 
343 
353 
36.1. 

363 

364 

365 

365 

366 

370 

371 
373 
400 
412 
424 
427 
434 
436 
441 
447 
457 

461 

462 
467 
505 
510 
522 

522 

522 


C 

C  STARTING  GUESS  IS  INCGOFRESSIBLE  FLOW 
00  50  1  =  1, NY 
DO  50  J  =  i , NZ MX 
50  RON ( I , J ) = 1. 0 

DO  26  NDsN=l,mEN 

CALL  MERMACH(NS,NR,NZ,H,T,NY,X,Y,  RON,  Z,  NCEII) 

DO  25  JS=1;NS 

J1=NZ (1 , JS)  SJ2=NZ (2 , JS I 

00  24  1=1, NY 

ITEST=0 

IFCI.NE.NYIGO  TO  1000 
1  =  1-1 
ITEST=1 

1000  HH=  H ( I, JS I 
TT=T(I, JS) 

D  Y= Y  ( I  ♦  1 )  -  Y (I) 

H0  =  HCI,ll-0.r- 
DO  23  J=J1,J2 
OZ=Z(J»l)-Z< J) 

IFCI.EQ.NY) 1=1-1 
JNYI=J*NY*I 

Kl= JNYI-NY  $  K2=K1*1  ?  <3=JNYI  IK4=K3*1 

X  1  =  X  ( K1 )  S  X2  =  X  ( K2)  Z  X3  =  XCK3>  7.  X4=X(K4I 

TlMl+C  (X4-X2)/0Z>  ♦•»2)/<X4*CX4-X3>  i  "2 

T2= ( 1 M IX4-X2)/0Z) •*£)/ (X2* CX2-X1I>"2 

T3=  1 1*(  CX3-XD/OZ)  **2)/  (X3*  CX*,-X3>  >  "2 

T-»=  ( 1+  {  CX3-X1I/IJZ)  "<.S/Ui“  (X2-X1))  "2 

TS=<Ti*T2*T3'T4)*0Y**2/4.U 

IFCITEST.NE.1IGO  TO  1001 

I  =  N  Y 

X1=X2 

1001  CONTINUE 
KR=  0 
QT=1.C 

OU  21  MNIT=1,MNMAX 
QTC= A  3S ( CT) 

EPSM=RON<I,J> 

A1=(2.'HH-TT/X1"2>*R0N(I,  J)  »*2 
A2=2.*HO*C'ONCI,J)  "  ( GA 'HI.  ) 

A  3= A2* ( GAM+1 . ) 

QT=C A1-A2-TS)/C2.*A1-A3» 
aT=OT*°F 

IF(OT».LT.A'JSCCT)  SGO  TO  211 
R0N(I,JI=R0!11I,JIM1.-0T> 

IFfADSIEPSH-RONtl, Jll.LE.ETOLIGO  TO  22 

21  CONTINUE 
GO  TO  Z 12 

211  TVAR=(1.0'TT/IX1"2) )/2,0 

RON  ( I ,  J  I  =5Qr'T  C  C  CHH-TVAR>/H0l  "5> 

RPM-’CN  C I , J} 

W  ° I T  E (6,204) I, J,PPH 

204  FORM4TC5X, 'DENSITY  ADJUSTED  TO  STABILIZE  COUf-LcO  PROBLEM*, 5X, 
1*1  =  *, 12, *J  =  *, 12, 'DENSITY  Sc.T* ,E12.4 1 

212  CONTINUE 

22  CONTINUE 
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522 

530 

541 

544 

546 


562 


567 

570 

575 

576 
602 
604 
61Q 
614 
620 
622 
630 
63  5 
647 
657 
666 
676 
706 
714 
736 
747 

751 

752 
764 
776 

1000 

1012 

1024 

1036 

105u 

1052 

1055 

1057 

1061 

1065 

1067 

1072 

1100 

1106 


111C 

1111 


C 

c 

c 

c 

c 

c 


23  CONTINUE 

ir(J‘J.EO.NS»  R0N(I,NZNX)=E0N(I,J2J 

24  CONTINUE 

25  CONTINUE 

CALL  MERI!ACH(NS»N9,NZ,H,T,NY,X,Y,R0N,Z»NDEN) 

***♦»»♦»**•*» '!»»»»*■*. **<:»•*»*«*»***«■  •*** 

BEGIN  THE  NEWTON-RAPHSON  ITERATION  * 


DO  20  IT-1» ITMX 

EVALUATE  THE  RESIDUAL  VECTOR  ... 

DO  2  K= 1 , NOD  E 
£(K1=C. 

SUM  THE  ELEMENT  CONTRIBUTIONS  TO  RESIDUALS 
DO  5  JS=1,NS 

J1=NZ(1,JSI  Z  J2=NZ  <2. JS) 

DO  4  1=1, NY1 
DY  =  YCI«-1)-Y(I) 

HH=H(I, JS) 

TT=T( I, JS» 

DO  3  J=J1,J2 
OZ=Z(JH>-Z(J> 

K1=JNYI-NY  i  K2=Ki*l 
X1=X(K1»  i  X2=X(K2) 

T 1= ( 1 ♦ ( (X4-X2J/DZ) **2)/(X4*«X4-X3»)  **2 
T2=  (1 K <X4-X21/DZ)**Z)/(X2V(X2-X1M **Z 
T3=(1M  (X3-X1>/0Z>**2)/(X3MX4-X3»I  **2 
T4=  lit!  (X3-Xl)/OZ>  **21/ (XI* IX2-X1) J  **2 
TS=( Tl+T2*T3 ♦Tt) *0Y*«2/4.Q 

RN=  (RON  (I ,  Jl  ♦RON(I,J+l)  VriON  1 1*1 ,  J  J  ♦  RON  < 1*1 ,  J»1 »  1/4.0 
IF ( J.EQ. J2I RN= (RON (I, J) ♦ROM <1*1, J»  I/2.C 

RN2=RN**2 

TS=TS/RN2 

SS=  (1.0/Xl*»2+1.0/X2**2fl.2/X3**2*l.L/X4**2t/4.0 
FF= (GAM-1. C J  *HH/GAN+(TSfTT*SS)/(2.0« GAM) 

PN=°E/RN 

E(K1)=E (K1)-0(X4,X3,X2*X1) 

E(K2)=t (K2)*U»X3,X4,X1,X2) 

E(K3)=E  (K3)-D(X2,X1,X4,X3» 

E  < K-* »  =E  (K4)  *D(X1,X2,X3,X4) 

CONTINUE 


JNY 1=  J*  NY* I 
K3= JNYI 
X3=X(K3» 


$  K4=K  3* 1 
5  X4=X(K4) 


CONTINUE 

CONTINUE 

MODIFY  RESIDUALS  AT  SPECIFIED  HALL  NODES 
IF(NF.EQ.O)  J?=NZMX 
IFtNF.EU.il  J2=N  Z ( 1 , NSI 
DO  6  J= 1 , J2 

JU=J*NY  I  JL= JU-KY+1 

E(JU)=X<JU»-G<J> 

E(JL)=X(JLI-FIJ) 

CONTINUE 

RESIDUAL  NORM  AND  CHECK  FOR  CONVERGENCE  ... 
EMAX=C. 

OO  7  K=1 , NODE 


RSOLVE 


47 


1 


RUN  VERSION  KE9  74  16147  L'6/24/74 


1112 
1123 
1126 
1132 
1132 
114  0 
114C 
1162 
1204 
1204 
1204 


1212 

1214 

1215 
1224 

1227 
1230 
1234 
123b 
1242 
1246 
1252 
1254 
1262 
1267 
1331 
1311 
1320 
1313 
134C 
1346 
1370 
1401 
1403 
14  J  4 
1416 
1430 
1432 
1450 
1466 
1504 
1522 
1542 
1561 
1576 
1614 
1633 
165  J 
1666 
1704 
1707 
1711 


7  EHAX= AHAX1 (t WAX, AOS <E <K )  )  ) 

IF<MOD(NOEN,NP&T).NE.O>GO  TO  777 
WRITE  <6 ,299) 

299  FO^MAT(lHi) 

WRIT*.  (6,3301  NO  EM 

300  Ff)»'MTI3(/»  ,5X, ’ITERATION*. 13, 5X, 'SOLUTION  AMO  EPRDP  P"SIDUAL*) 
W p 1 1 L  <6,3 31)  (X  ( K ) «  K  =  1 » NODE) 

WRIT!-:  (6,301 »  <;;<K>,  K=l,  NODE) 

301  FORMAT  < / , ( 5 X , H E 1 1 • 3 ) ) 

777  CONTINUE 

IF(CMAX.LE.E70L)GO  TC  31 
C 

C  EVALUATE  THE  JACG3IAA  ELEMENTS  NUMERICALLY  ... 

00  9  K=  1, NOD E 
OO  A  L=1,N3AN 
B  C  <  K  ,  1. 1  =  0  • 

9  CONTINUE 

C  SUM  THE  ELEMENT  CONTRIBUTIONS  TO  JACO0IAU 

00  1?  JS= 1 »  NS 
J1=NZ(1,JS)  I  J2=NZ (2, JS> 

00  11  1=1, NV1 
OY=Y(If l)-Y(I) 

HH  =  i<  C I, JS) 

TT=T ( I, JS) 

00  1U  J=J1, J2 

OZ=2(J*-l)-ZIJ)  T.  JNYI=J*NY«-I 

K1  =  JNYI-AY  Z  K2  =  K1U  S  K3  =  JNYI  £  K4  =  K'3*1 

X 1  =  X ( K1 )  :  X 2  =  X  ( K2 )  >  X3=X(K?)  £  X4  =  X(K4) 

T1  =  ( 1  ♦  (  (X4-X2)  /0  7. )  **2) / ( X4*  IX4-X31!  tf¥2 
T2=  <1 «■  C (X4-X2) /DZ) **2)/ ( X2* (X2-X1 ) ) **Z 
T3=(1M  (X3-XD/OZ)  **2)/<X3*(X4-X3  ))‘!'' 2 
T4=  ( 1  ♦■  (  (X3-XD/DZ)  ** 2*  /  (XI*  (X2-X1)  )  **Z 
TS=  (T1*T2*T3 *T4) *0Y**2/4.3 

RN=  (RON  1 1 ,  J)  *EON  ( I.JUI  f  RON  ( I U  ,  J)  *  RON  (I*  1 ,  J*1 )  >/4.0 
IF(  J.F.Q.J2)RM=(PON(I,J)  »RGN  ( I  ♦  1 ,  J  >  >  U  .  C. 

RN2=PN**2 

TS=TS/RN? 

SS=  (l.l/Xl**  2H.C/X2**2n.  :/X3**2*l.i,/X4**2 1/4.0 
FR=(GAM*l.(i)  *HH/GAMf(TS*-TT*SS)/  (2.ri*GAM> 

PN=  PE/Q  N 

013=D(X4,X3,X2,X1>  %  02u=0 ( X3 , X4 , XI , X2) 

030=0 <X2, XI, X4,X3)  £  U 4C =P < X 1 , X 2 , X3 , X4 > 

C(Kl,l)=C(Kl,l)-(0(X4,X3fX2,Xl*0X)-Jir,  )/OX 
C(Xl,2)=C(Kl,?)“(OlX4,X3,X2«-nX,Xl)-Ul!.(/DX 
C  <K1  ,NYf  1)  =C  (K1,NY  U)  -  (0  (X4,  X3*0X,X2,M)  -012)  70X 
C(K1,NY*2)=C(K1,NY*2)-(0(X4*-0X,X3,X2,X1)-01C)/UX 
C  ('<2, 1)  =  C(K2 ,1)  ♦  (0  (X3,X4,X1 ,  X2  *  Q  X  >  -  f) 2  3  )  /  D  X 
C(K2,NY)  =  C(K?,MY  )  (0(X3  +  r>X,  X4  ,  X  l ,  X2)-Q 231/OX 

C(K2,KY*1)=C(K2,NY*1!  MOIX3  ,X4»DX  ,X1  ,X2)  -020 /OX 
C  (K3,i)=C(K3,l) - (0 (X2,Xi, X4,X3*0X)-0cG ) /OX 
C(KT,2I  =C(K3,2  )-(0  (X,?,  XI,  X4*0X,  X3)  -03C  )  /  OX 
C(K4,  1)=C(K4,1) «-  (0  (X1,X2,X3  ,X4»OX)-04C  )/0X 

10  CONTINUE 

11  CONTINUE 

12  CONTINUE 

C  MODIFY  JACODIAN  AT  SPECIFIEO  WALL  NODES 


RSOLVE 


48 


RUN  VERSION  FE9  74 


16*47  66/24/74 


1714 

IFINF.EQ.O)  J2-N  Z  HX 

1716 

IF(NF.EQ.l)  J2=N7(i,FS) 

1722 

00  14  J=1,J2 

1724 

JU=J*NY  %  JL- JU-NY^l 

1727 

C(JU,.)  =C(JLtl>=l. 

1737 

00  13  L  =  2»  N9AN 

1740 

C(JU,L)=CUL,L)=0. 

1750 

13 

CONTINUE 

1753 

C  ( JU- 1 »  2) =3 • 

1757 

IF( (NF. lC.1I  .AMO. 1  „.£0. J2)  )  GO  TO  14 

1766 

C ( JL 1 1 1  NY)  =  0 • 

1771 

C(JU-1,NY*2I=C. 

1775 

p 

14 

CONTINUE 

u 

c 

SOLVE  THE  SYMMETRIC,  POSIT I ','E-OEF  INI!  r  JACOBIAN  SYSTEM 

c 

USING  A  BANOtO-n OLESKY  BE COMPOS IT  ION  ... 

2C0C 

c 

CALL  BANSOL INBAN, NODE, NR, C,E> 

c 

c 

NEW  ESTIMATE  OF  THE  SOLUTION  VECTOR  ... 

2004 

OXMX=C. 

2005 

00  15  K=i , NODE 

2012 

X(K)=X(KI-E(K> 

2016 

OXMX  =  A,!AXl(OXHX,  A3S<£IK»> ) 

2025 

p 

15 

CONTINUE 

2027 

p 

20 

CONTINUE 

2032 

26 

CONTINUE 

20  34 

WRITE (o ,232) 

2o4C 

31 

CONTINUE 

20  4  0 

NDEN=Ii)EN 

2042 

CALL  ME RK  ACH  *NS  ,  NR,NZ  »H,  T  ,  ,‘J  Y  ,  X  » Y  ,  RON ,  7  ,  NDEN) 

2061 

201 

FORMAT)*  I HP ,  OFF 0  ARFAY  SIZE  FOR  EXECUTION  *, 

*  NR, NS, NY,  f.ZNX.NOOE  *,515) 

2061 

p 

2C2 

FORMAT)//, 15X,*SCLUTI0N  OIO  NOT  CONVERGE*) 

2061 

RETURN 

2C62 

END 

RSOLVE 
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15 

15 


17 

22 

24 

26 

27 

40 

45 

45 

72 

76 

100 

102 

104 

106 

107 

137 

143 

143 


SUBROUTINE  GUESS (X ,F  ,G , NZMX  ,NY , Y , XT , XH , Y T , YH) 
OIHENSI ON  X(NY,1» ,Yll> ,F(1» ,G(1) 

NY1=NY-1 
C  ENO  R  VALUES 

C  PRESCRIBE  SURFACE  PS  I  VALUES 
YH=-<XH**2»/2. 

YT=-<XT**2»/2. 

C  CALCULATE  INTERMEDIATE  PSI  VALUES 
V(1)=YH 
DO  10  1=1, NY1 

R  =  XH*  (XT-XHJ  *FLO AT  (  D/FLOAT  (NYl) 

10  Y<m»=-0.5*R*R 
PETURN 
ENTRY  COORD 
XD=<XT/YH1**2-1. 

DO  15  J= 1 , NZ  MX 
A1=F{ J) ♦*? 

A2=G ( J) *  *2 
91= A2- A 1 
00  15  1=1, NY 

X(I,J)  =  F<J>MSGRT<-2.*Ym>-XH>MG<  J) -F  (  J  >»  /  I X  T -XH) 
15  CONTINUE 
RETURN 
ENO 


GUESS 


50 
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3 


3 

3 

17 

22 

23 

37 

61 

43 

44 

45 
45 


45 

65 

65 

105 

105 

125 

125 

144 

144 

146 

147 
152 


153 

165 

173 

173 

231 

201 

207 

207 

216 

22C 

237 

237 

242 

243 
26  3 


263 

275 


PROGRAM  TC0MFL0IINFUT,0UTPUT,TAPE5*INPUT,TAPE6=0UTPUT) 

OIMENSION  NZI2.3J ,FI24) ,GI24) ,H 16,31  ,T (8,3) 

1, XI 8, 24 1 ,YI8> ,ZC 24), Cl 192,141  ,£(192) 

1, RON I  9, 24) 

DIMENSION  ALFA  13  )  ,  EcTAIi)  , OMEGA  13) 

666  REAO(5,101INV,NZMX,NS,NF 
IF  I  EOF, 5) 1,2 

1  CALL  EXIT 

2  WRITF(6,1031INY,NZXX,NS,NF 
PI=3. 1415927 

NY1=NY-1 
NOOE=NY*NZMX 
NP.=  192 

101  FORMAT (415) 

1001  F0RMATI1H1,3I/) ,15X,*NUM9r?  OF  STREAMLINES  =*,I2, 

1  /,15X,*NUM3ER.  OF  GPIU  3i.GMTS  =*,I2, 

1  / ,  15  X , *NUMGER  OF  FLOW  SECTNS  =*,I2, 

1  /,15X,*3=ENCL0SED,  1=FRE«.  EXIT*, 12) 

REAO (5,  1021  IJ.ZIJ)  ,F  IJ) ,G(J) »J=1»NZMX) 

102  FOPMATII5,3F10.5) 

HRITEI6, 1002) IJ.ZIJ), FIJ), GIJ) ,J=1,NZMXI 
1302  F0RMATI2I/) ,5X,*S£ GHENT*  »5.<»* ORDINATE*  ,5 X,*HU3  SURF ACE*,  5 X, *T IP  SU 
1RFACE*,  /, 15 X, I4,6X,F13.5,5X,F1C.5,5X,F10. 5)  > 

PEAO  15 , 137)  IUZIl,J),NZI2,j) ,J  =  1,MS) 

107  FORMAT! 3121101) 

WRITE  16, 1037) IINZ 1 1 , J) , 1=1 , 2) , J=1 ,NS) 

1007  FOR'  AT  13  I/) ,5X, ‘SECTION  DELIMITERS* ./ , I 3X ,2 IIS ) J 
C  PRESCRIBE  SURFACE  PSI  VALUES 
XH=FI1) 

XT=G I 1 ) 

YH=-3.5*XH»*2 

YT=-„.5*XT**2 

C  CALCULATE  INTERMEDIATE  PSI  VALUES 
C  SET  HU 9  AND  Til*  COORDINATES 
C  SET  PSI  LINE  VALUFS 

CALL  GUESS  I X ,F ,G , N ZMX.NY , Y , XT , XH, YT , YH) 

C  ENTHALPY  ANO  ANG  MOM  TERMS 
REAOI5.1P3) IFLAG 

103  FORMAT  1 1?) 

RE AO  1 5 , 1 341 1 HCONS 
1041  FORMAT  I Fi 0«  5) 

WRITr  16 , 1042  I HCONS 

lli42  FORMAT!//, 5X,*INCTCENT  MACH  NO.  *,F5.2) 

HCONS= 15.3/ I2.0*HCCMS**2) ) *11. C*3,2*HCONS**2) 

777  IF! IFLAG. EQ.11GO  TO  888 

REA0I5, 1 J4) I  ALFA! J)  ,  BET  A  I J) .OMEGAIJI ,J=1,NS) 

1C4  F ORM A T I  3  c  13 . 5) 

IF  (EOF, 5)3, 4 

3  GO  TO  066 

4  WRITE  16, 1003)  IJ.ALF  A  IJ)  ,BET  A  I  J)  ,  OMEGA  (  J)  ,  J=1,NS) 

1033  cORMAT  I  3  (/  I  ,  SX  , »  SECT  ION*  ,r>X  ,  *  ALF  A-C  O..F  *  ,  5X  ,  *6E  T  A- COEF*  ,5X, 

1* OMEGA  *,/,l5X,I5,3l5X,Fl 3.5) )) 

C  STARTING  GUFSS  IN  X 

CALL  COOR O I  X , F , G , NZMX , NY , Y , XT , XH, VT , YH > 

OO  35  1=1, NY1 


TCOMFLO 


.7 


RUN  VERSION  FEB  74  16147  06/24/74 


277 

H  (  I  ,  1  )  = HCONS 

361 

T  ( 1, 1 )  =  0.C 

30  2 

00  35  J=2,NS 

304 

T  II  ,J>  =  (6.5*  ALFA(J)*  (Y(I+1>«Y(  I)  >  -9£TA( J)  )  **2 

314 

OR  V  =  SaRTlT(I,J>>-SCRT(TCI,J-l>) 

326 

35 

H ( I f  J )  =  H  < I ,  J- 1 )  ♦  CME0A(J)*U RV 

341 

GO  TO  41 

34  1 

888 

RF AO (5, 105) C OMEGA (J)  ,J=1,NS) 

364 

IF (-OF,  5)6,7 

357 

10  5 

FORMAT ( 8P10.5I 

357 

6 

GO  TO  666 

360 

7 

RE AD ( 5 , 10  6)  ( (T(I,J) ,  1=1, NY)  ,J=1,NS) 

377 

106 

F  OR  MA  T ( 8F 10,5) 

377 

WRITE (6,10.6)  CJ, ONEGA ( J) , (T (I,J) » 1=1,  NY)  ,J=1,NS» 

422 

13  06 

FORMAT  (3  (/) ,5X,*SECII0N*,5X,*ANG. VEL.*  ,5X  ,  *  MOMENTUM*  ,/ , 
1(5X,I4,3X,F13.  5,3X,ri2,5) ) 

422 

00  36  1=1, NV1 

4?': 

H  ( I  ,  11 =HCONS 

426 

DO  36  J  =  2, NS 

427 

ORV  =  SQRTd  (I,J))-SQRT(T(I,J-1)> 

441 

36 

H ( I , J )  =  H(I,J-1)  ♦  OMEGA ( J) *ORV 

453 

41 

CALL  SECOND (Tl) 

455 

CALL  R3CL/E(NF,NP,r.S,Ny,NZ,C,E,F,G,H,T,PE,X,Y,Z,FON) 

475 

CALL  StCONO(T2) 

477 

WPITE (6,203) 

503 

203 

FOP MAT ( 1H1,5 (/)  , 23 X, ’SOLUTION  POINTWISt  BY  VERTICAL  SECTIONS*) 

503 

00  4C  J= 1 , NZ MX 

605 

WPITE(fc,20i) J , 7 ( J) 

514 

WPITE(6, 205) 

520 

2C5 

FORMAT  (/  ,  14X,*LINE',3X,*STREAMFUt!CT10N*,5X,*P6DIUS*  ,  5  X,*  AXIAL  VELO 
1C.* ,3X, *T ANG.  VELOC.*) 

520 

DO  40  1=1, NY 

522 

IF(  IFLAG.EQ. 1) GO  TC  37 

524 

DO  33  M  =  1 , NS 

525 

IF(J,NE.NZ(1,M))G0  TC  33 

530 

K=M 

531 

33 

CONTINUE 

534 

V  =  -( ALF A( K) *Y(  I) *B£T  A ( K ) )/ (X(I,J)/XH) 

544 

GO  TO  371 

545 

37 

V=-S0RT ( T ( I , K ) ) 

552 

371 

IF ( I,  NL  .  NY) GO  TO  3  8 

554 

P=Y (MY) -Y (N v -1 ) 

556 

Q=Y (NY) -Y (NY-2) 

560 

DE=P*  Q* CP-O) 

563 

A1=P**2 

564 

A2=Q**2 

565 

c\» 

<1 

• 

(l 

r 

<r 

567 

Fl=X(NY-2, J) **  2 

572 

F2=X(NY-1, J) **2 

574 

F  3=  X ( NY , J)**2 

576 

W  =  -2.*OE/  (-iU*Fl*A2*F2»A3*F3) 

605 

H=W/ROH ( I, J) 

611 

GO  TO  39 

611 

38 

CONTINUE 

611 

W=-2#MV(I*1I-Y(I)  )/(X(I*l,J)**2-X(I,J)**2) 

623 

W=W/RON (I, J) 

TCONFLO 


52 


RUN  VERSION  FE  3  74 
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626  39  WRITE(6,2G2>  I,YfI)  ,X(I,J),W,V 

646  40  CONTINUE 

653  201  FORMAT (//  ,2X.*5EGM£NT*«I3|5X»*0R0INAT£  Z(J)*,E12.4> 

653  202  FOR‘1AT<13X,I5,5X,tl2.4,4X,El2.4,2X,E12.4,3X,tl2.4l 

653  TIM=T2-T1 

655  WRITE(6»?j4)TlM 

663  204  FORMAT < 3  </>  ,6X,*S0LN.  TIM£*,F7. 3,*  C.P.  SECONDS*) 

663  GO  TO  777 

664  END 


TCOXFLO 
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16 

16 

16 

16 

20 

21 

25 

27 

33 

36 

40 

41 
43 
45 
53 
65 
71 

105 
112 
121 
122 
122 
123 
126 
134 
136 
14  C 
142 
145 
145 
145 
150 
152 
163 


163 

171 

176 

176 


S'JBROJT INc  MERMACH <NS,NR,NZ,H,T ,NY,X, Y, RON,Z,NDENI 
DIMENSION  NZC2.1I ,H(NY,1I  ,T (NY, 1) ,XC1>  ,  Y(l> ,Z( 1) , RON t NY ,11 
DIMENSION  \IP  (6,42) 

REAL  MMN , MSM 

C  MERIDIONAL  MACH  NUMBER  CALCULATION  BEGINS 
NY1=NY-1 
OO  43  JS=1,NS 
Jl=  NZ  (1  ,  JSI  iJ2  =  NZ(2,JS) 

00  43  1=1, NY1 
HH=H( I, JS) 

TT=T ( I, JS) 

K=0  $  L=1 

VMMN=G. u 

00  39  J=J1,J2 
JNY I  =  J*  NY  +1 

Kl= JNY I-NY  SK2=K1H  £  <3=JNYI  2K4=K3*i 
X 1=  X ( K1  I  IX2=X<K2)  £X  3=  X ( K3 )  «X4=X(K4I 

V=SQR1 ( TTI/Xl 

H=-2.  *  C  Yf  IU)-Y(I)  )  /  (X2**2-X1**2I 

H=W/ EON  (  I,  J) 

u=-w*(X3-xii/(Z(j«-ii-Z(jn 
U  =  U**2 
H  =  W  *  *  2 
V=V**  2 
TOT  =  U*V  *-H 

MMN=2.5*  <UfWI/'HH-3.5*TOT> 

VP< I, J> -MMN 

MMN  =  AoS C  MMN I 

IF( VMMN.GE.MMNIGO  TO  38 

K=I  ?  L=J 

VMMN-HMN 

38  CONTINUE 

IF( MMN. LE.1.IG0  TO  39 
MEM=S0RT ( MMN) 

WRITS (6 , 1C  2) I, J,MEM 

102  FOpM  A  T ( 5X,*tNCCUNTC?tO  MERIDIONAL  MACH  NUMBER  IN  EXCESS  OF  UNITY*, 
15X , *  I- * , I2,*J=*,I2 **  MER.  MACH  NO.  =*,F5.2) 

39  CONTINUE 
4C  CONTINUE 

C  MERIDIONAL  MACH  NUMBER  CALCULATION  ENOS 
RETURN 
END 


MERMACH 
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C 

C 

C 

C 

C 

C 

c 

c 

10 

c 

c 

1C 

11 

15 

20 

22 

24 

30 

35 

46 

5u 

54 

56 

6C 

64 

7C 

72 

74 

111 

121 

C 

c 

124 

125 
127 
132 
134 
136 
147 
156 
160 
162 
164 
167 
171 
173 
20  o 

C 

216 

216 


SUBROUTINE  0  ANSOL  (  H ,  N  ,NR  ,  A  ,9) 

SOLVES  A  3ANCF0,  SYMMETRIC*  POSIT IVE-OEEINITE  SYSTFM  CT 
N  LINEAR  EQUATIONS  USING  PACKED  S 1  OR  AGE  AND  A  CHOLESKY 
DECOMPOSITION  FOLLOWED  3Y  EORWA&n  AND  HACK  SUBSTITUTION 
NR  IS  THE  ROW  OIMENSION  OF  A(I,J>  WHEN  IT  WAS  STORlD. 

H  IS  THE  HALF-°ANQ  SIZE  (SUPERDIAGONAL)  PLUS  ONE. 

ON  ENTRY  THE  SINGLE  c IGHT-HAND  VECTOR  IS  STOR'D  IN  B(D 
AND  O(I)  CONTAINS  THE  SOLUTION  VrCTOP  ON  EXIT. 

DIMENSION  A ( NR  *  1 ) *  3 ( 1 ) 

CHOLESKY  DECCMPOSI TION  WITH  PACKED  STORAGE  ... 

OO  5  1=1, M 
SUM= A (1,1) 

KMX=MIN0 (I,M) 

IF ( KMX. LT.2 )  GO  TO  2 
00  1  K=  2  ,  K M X 
TMP=A(I-K»1,K> 

1  SUM=SUM-TMP*TMP 

2  A (1 ,1 )  =  RA=SORT (SUM) 

RA=1./RA 

JMX=;1INC  (H,H-IH) 

IFIJMX.LT. 2)  GO  TO  5 
DO  4  J=  2  ,  JM  X 
SUM=  A ( I , J) 

KMX=MINO(I,M-Jfl) 

IF ( KMX. LT.2 )  GO  TO  4 
DU  .3  K=  2  ,  KMX 

3  SUM  =  SUM-A(I-KH,K)*A(I-K+1,  J*K-1) 

4  All, J)=PA*SUM 

5  CONTINUE 

FORWAFO  AND  BACK  SUBSTITUTION  WITH  PACKED  STORAGE  ... 

00  7  1=1, N 
SUM=  3 ( I  ) 

K  M  X  =  M I N  0  (I,H) 

IFiKMX.LT, 2)  GO  TO  7 
00  6  K=  2, KMX 

6  SUM=SUM-A(I-K*1,K)*B«I-K»1) 

7  3 ( I ) =  SU  M/ A  1 1  , 1 ) 

DO  9  11  =  1, tl 

I  =  N“ 1 1 ♦  1 
SUM=B(I ) 

KMX=MIMCC1I,M> 

IF (KMX. LT.2)  GO  TO  9 
DO  2  K=  2 , KMX 

8  SU“=SUM-A(i ,K) «B(I*K-1) 

9  3 ( I ) =  SUM/ A ( I , 1 ) 

RETURN 

END 


8 ANrOL 


55 


NUMBER 

OF 

streamlines  - 

8 

NUMBER 

OF 

GRID 

SEGMTS  = 

24 

N'JMC.ER 

OF 

FLOW 

SECT  NS  = 

3 

C=FNCLOSED 

»  1^ 

FREE  EXIT 

0 

SEGMENT 

OROINATE 

HUB  SURFACE 

1  IP  SUPFACl 

1 

-9.420C3 

1.00330 

3 .L  GCuv 

2 

-7.04G09 

1.0.-300 

3.0GGoP 

3 

-5. 13G  30 

1.00030 

3 . 0  r  c  c  r 

4 

-3.6CCC  3 

1.00033 

3.0CCCC 

5 

-2.38000 

1.CC200 

.3 . 0  u  0  C  * 

6 

-1.4  llJCO 

1.0CC9C 

3.G90GC 

7 

-.630  CO 

l.GOOOO 

3.0 0170 

8 

-G.OPGOO 

1.03033 

3  .  G  0  65* 

9 

•  500  Cl 

1.01009 

3. G 1720 

10 

1.00000 

l.tUOQO 

3  .  G  2  6  3  (. 

11 

1.5CG  tO 

1.00C03 

3.C410C 

12 

2.C0GQ0 

1.GC390 

3.05460 

13 

2.5CC 00 

l.OCJOO 

3.06810 

14 

3 . 0  0  0  0  j 

I.GCOjO 

?.09cer 

15 

3 . 5C  0  0  0 

1.CC030 

3.G919C 

16 

4.00000 

1.0CC33 

3.  IC06L 

17 

4 • 630  CO 

l.CCOJJ 

3. 1074* 

18 

5.41000 

1.C0J00 

3.1086C 

19 

6 . 3  8  u  C  0 

l.onojo 

3.1086C 

20 

7.600  DC 

1 . 0  u  C  «*'  0 

3.10860 

21 

9. 130  CO 

1.0  CO  JO 

3.1086C 

22 

11. 0  40  CO 

1 .  CIO  3  9 

3.10660 

23 

13.42000 

1.00000 

3.1086C 

24 

16.4COvO 

1.00000 

3.1086C 

SECTION  0 

"LIMITERS 

1 

8 

9 

16 

17 

23 

INCIDENT 

MACH  NO.  .20 

SECTION 

ALFA-CGEF 

BET  A-COEF 

OMEGA 

1 

O.GQCOU 

O.CCJPO 

C  .GOCOC 

2 

C. DLOGO 

2.G00G3 

D.CCuOO 

3 

O.GGCGO 

C. 00300 

C  .OOGOC 

SOLUTION  FOINTWISc  VEPTICAL  SECTIONS 


segment 


SEGMENT 


SEGMENT 


SEGMENT 


SEGMENT  5 


OROINATE  Z<J)  -9.42CG£*C0 


LINE 

STREAMFUNCTION 

RADIUS 

AXIAL  VELOC. 

T  AUG.  VELOC 

1 

-5.CC00E-91 

1.  G’.CE  +  OQ 

9. 9735E-31 

“0. 

2 

-8, 2653E-U1 

1.28GGi>G3 

9.9fcfi9E-0l 

-9  . 

3 

-1.2347E«-Ci: 

1.572oE»u9 

9. 9982E-G1 

-0. 

A 

-1.724F£M'j 

1.  8533/1  *-00 

1.  Plo5E*ui. 

-p . 

5 

-2.2959EU’  G 

2  •  1 A  ‘V  r.  E  *■  C  P 

5 . 0 jC8ErQP 

-c . 

6 

-2.  9  4  9 1 E  ♦  f!  C 

2.  A293EUC 

1. Co  1 ?E 

-0  . 

7 

-  3  ,  6  8  3  7  l  ♦ .)  C 

2. ?1aoE*G3 

1  •  3  l’  1 3E  +  0  r 

-p. 

8 

-A.5G0CE*0i) 

3»u,/3iE*’uC 

1.0ul3E+Cu 

..r  , 

ORDINATE  ?<J)  -7. C4CQE+CC 


LINE 

STREAMFUNCTION 

RADIUS 

AXIAL  VELOC. 

T  ANG, 

1 

-5.3  ICCE-C1 

1 .  C  3  J  0  E  ♦  3  0 

9,973  At -Cl 

w  n 

2 

-8. 2  653  E-3 1 

1.2856E+30 

5 , 96  88E- C 1 

*  W 

3 

- 1. 2347E  «-G  3 

1 • 5726E*  53 

9. 5931E-31 

-0 

A 

-1.72A5EOU 

i»8533E*3J 

1.3GS5E+00 

5 

-2. 2959E*0C 

2. 143  8F ♦ 00 

1. COG  8E»Cw 

.0 

6 

-2.949PE«"JC 

2.4293E+G0 

1,331 3E ♦ G  C 

-c 

7 

-3.6837E*0C 

2.714oE*-C0 

1.33 13£ +CC 

-■ 

8 

-A, 5QCOE+OC 

3,  S3 0 Q E ♦ 03 

1.0G13E+0C 

-p 

ORDINATE  Z(J>  -5.13COE«-CO 


LINE 

STREAMFUNCTION 

RADIUS 

AXIAL  Vt'LOC. 

TANG.  VELOC 

1 

-5.GGCCE-C1 

l.C'JOE  +  GC 

9. 9735E-C1 

-C  . 

2 

-8.2653E-31 

1,  86GE  +  33 

9.9887E-01 

-0. 

3 

-1. 23A7E*  GC 

1.  j 726C  + 33 

9.  99 82t -G  1 

-3  . 

A 

-1.7245F* 30 

1, 8533t  +0C 

1.3335E*GP 

•  n 

5 

-  2. 2  959E  *o  u 

2.1A38E+3J 

l*QSj8E»00 

-3. 

6 

-2.9A90c«G3 

2. A293L+03 

1.03 13E+t3 

-3. 

7 

-3.6837EH1C 

2.7146G*03 

1 . Co 13E*G  C 

-1 . 

8 

-A.5COCC+OG 

3.0JG0c»C3 

1. JC  13E  +  CG 

-J. 

OROINATE  Z(J»  -3.6G00E»Q!) 


LINE 

STREAMFUNCTION 

PADIU3 

AXIAL  VELOC. 

TANG,  VcLOC 

1 

-5. JOOGt-Jl 

l.OGjjt^GQ 

9. 979AE-01 

-0. 

2 

-8. 2653E-3 1 

1. 2866E ♦ GO 

9. 9697E-01 

-7  . 

3 

-1.23A7E*0C 

1. 5726£ful 

9. 9992E-C J 

-G. 

A 

-i. 7  245C  *0  o 

1.85B2L+C0 

1  •  3 1!  35E  ♦ 'J 

-C  . 

5 

-2.2  959E* 3  3 

2.1438E+3U 

1.3309£*-0C 

-0  . 

6 

-2. 9A9n£»Q0 

2. 4292E«-uG 

1. 0013E+0G 

-n. 

7 

-3.  6  337E*  Q  3 

2*  7146E  ♦  00 

1.3C12E»cr 

-p . 

8 

-4.5C00E*0G 

3.GGJuE  +  Cii 

1 • Go 11E+G0 

-o . 

OROINATE  Z(JI  -2.38002»GC 

LINE  STREAMFUNCTION  RADIUS  AXIAL  VlLOC.  TANG,  VELOC. 
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SEGHcNT 


SEGMENT 


SEGMENT 


SEGMENT 


SEGMENT 


3 

-1.2347F»G0 

1.6291E*0u 

9.63J7E-C1 

-1.2277c*CC 

4 

-1.7245F»00 

1.9238E+0E 

9.63376-01 

-l.CT96£*GU 

5 

-2.2959c»ac 

2.  2164E*  jc 

9.6137E-G1 

-9.  0238F.-C1 

6 

-2.949!.fc*00 

2»5«fl 26 * OC 

9. 5983c* . 1 

-7.97376-01 

7 

-3.6e?7t»iia 

2.7999E*.3 

9.57G4R-C1 

-7.1432E-M 

8 

-4.536CEOQ 

3.3919E*CC 

9.5509E-01 

—6.  46856—01 

ORDINATE  ZIJI  4.00006*00 

LINE 

ST  RcANFUNCT ION 

RADIUS 

AXIAL  WELOC. 

TANG.  VELOC. 

1 

-5.  G .GGE- j 1 

l.n000E*33 

9.29936-01 

*2  •  r  V  JC  E*  CC 

2 

-8.2C536-01 

1. 32706*03 

9. 56966-01 

-1. 55 72- *29 

3 

-l.2347E*03 

1.62836*0:' 

9.5758E-01 

-1.2279->oC 

4 

-1.72456*03 

1.9250E+G3 

9.5615051 

-1.039?E*„0 

5 

•2.  29596  *CC 

2.2193E*CC 

9.534G-: -61 

*9. 312l£-1 1 

6 

-2.9-*9:£*GC 

2. 513£E*uG 

9.51296-01 

-7.9586F-L1 

7 

-3.6C37F*10 

2. 9 Jb5£*0L 

9.47440C1 

-7.12636-Ci 

8 

•4.  5  CCCf  *0  0 

3. lb  066*03 

9.44976-51 

-6.4  5  54  L-  Cl 

ORDINATE  Z(J>  4.6303d 

*00 

LINE 

ST  REAMFUMCT ION 

RADIUS 

AXIAL  VELCC. 

TANG.  VELOC. 

1 

-5.0JiiCF.-31 

i.coooe*u: 

8.7992c-ul 

-C  • 

2 

-8. 26536-01 

1.  31 93  c.  *00 

9. 11926-El 

-0. 

3 

-1.2247EU6 

1.62256* to 

9.22156-ul 

4 

-1.  7245c  <  Cl  C 

1. 9215c*  G  ." 

9.27j9c-tl 

-0. 

5 

-2.29556*00 

2. 2104E+0C 

9. 2889" -01 

-a. 

6 

-2.949tE*0Q 

2.5147c*.! 

9.29666-01 

-0. 

7 

-3.6'337F*C0 

2. 3107t*C0 

9.27426-0.1 

_  *> 

•  • 

8 

-4.5CCCt*30 

3.IC74£*uC 

9. 2629c -C  1 

-c. 

OROIHATE  ZCJI  5.4133E*0Q 

LINE 

STRcAMFUNCT ION 

RADIUS 

AXIAL  VELCC. 

TANG.  VELOC. 

1 

-5.;:.!jrE-2l 

l.C0JQc*U0 

9. 0714E-C1 

_f» 

u  • 

2 

-8.  265.3E-3  1 

1 . 31 G  3c*  0  J 

9. 1 497 E- Cl 

-0. 

3 

-1.23476*00 

l.Gl*9E*30 

9.1936E-&1 

-0. 

4 

-1.72456*03 

1 . 9159t*jC 

9. 22226 -fc  1 

-c. 

5 

-2. 2 95 9c* CO 

2«2151E*QC 

9. 2363c- t 1 

-0. 

6 

-2. 949;.t  *  00 

2. 51336*0. 

9. 24656-t 1 

-0. 

7 

-3. b637t  *u  G 

2, 8109c*  GC 

9,2399c-!'! 

-3. 

6 

-4.5lia:F*03 

3.  13  36  c* . 9 

9. 2366E-C1 
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SECTION  VI 


EXAMPLE  RESULTS 

1.  CALCULATION  OF  DESIRED  AREA  CHANGE 


In  the  first  series  of  examples  considered  below,  a  stator  rotor  pair  are 
considered  to  exist  in  an  annulus  in  which  the  area  is  expanded  to  make  the  out¬ 
let  "average"  Mach  number  approximately  equal  to  the  inlet  average  Mach  number. 
Assuming  for  the  moment  that  a  single  Mach  number  exists  far  downstream  (as 
would  be  the  case  for  free  vortex  machinery)  it  follows  from  consideration  of 
continuity  for  isentropic  flow  that 
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or,  in  the  case  where  M,  =  M  , 
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Here  we  calculate  the  desired  area  ratio  using  this  formula,  but  consid¬ 
ering  the  enthalpy  to  be  a  mass  averaged  enthalpy  such  that 
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For  convenience  of  interpretation  the  examples  considered  herein  were  all 
^  iken  to  be  of  the  linear  form  considered  in  Ref.  3.  That  is  to  say  the  angular 
ntum  in  all  sections  is  of  the  form  YV  =  -Jt'f  +  8.  Thus  it  follows  in  this 
case  of  a  stator  followed  by  a  rotor  that  removes  the  swirl,  that: 
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area  ratio  is  the.i  calculated  from  Eq  (56). 
that  the  hub  radius  remains  at  unity  and  the 


It  is  assumed  in  these 
tip  radius  in  the 
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expanding  section  of  length  L  is  given  by 
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In  the  examples  considered  in  this  report,  L  was  taken  equal  to  5.4063. 


2.  FREE  VORTEX  FLOWS 


Figs. 16  through  21  show  the  results  of  the  calculations  for  flow  through  a 
stator  rotor  combination  in  which  a  free  vortex  swirl  is  introduced  by  the 
stator  and  removed  by  the  rotor.  In  all  cases  the  outlet  area  is  adjusted  to 
compensate  for  the  reduced  stagnation  enthalpy  as  described  in  the  previous 
section.  The  figures  are  arranged  in  order  t'  illustrate  the  effect  of  in¬ 
creasing  inlet  Mach  number,  other  pertinent  parameters  being  supplied  on  the 
figures . 

It  is  readily  apparent  that  the  perturbations  in  axial  velocity  brought 
about  by  the  effect  of  compressibility  increase  substantially  over  the  range  of 
Mach  numbers  investigated.  (For  example,  just  following  the  rotor  =  0.95 

for  Mq  =  0.2,  whereas  =  0.68  for  MQ  =  0.5).  The  overall  effect  of  the  reduc¬ 
tion  in  axial  velocity  at  exit  due  to  the  reduced  stagnation  enthalpy  is  also 
evident  in  the  same  figures.  (W  =  .98  for  M  =  0.2,  W  =  .91  for  M  =  0.5). 

It  is  obvious  that  in  these  cases,  where  the  vorticity  is  zero,  that  the  outlet 
velocity  must  be  uniform,  and  it  is  a  simple  check  on  the  accuracy  of  the  pro¬ 
gram  to  observe  the  predicted  uniformity  of  the  outlet  velocity  profile. 

It  is  interesting  to  note  that  in  the  calculation  of  these  free  vortex 
examples,  no  spurious  values  of  meridional  Mach  Numbers  in  excess  of  unity  ap¬ 
peared  throughout  the  iterations  in  the  calculational  procedures.  This  was  true 
even  though  in  the  case  where  MQ  =  0.5,  meridional  Mach  number  on  the  hub  just 
following  the  stator  was  M  =  0.74.  (The  total  Mach  number  at  this  location  was 
1.36).  In  the  examples  to  follow,  the  very  important  effect  of  the  vorticity 
in  distorting  the  axial  velocity  profiles  becomes  evident.  The  effect  of  the 
distortion  is  to  actually  introduce  high  Mach  numbers,  as  well  as  to  make  the 
eventual  streamline  pattern  more  distant  from  the  original  guess. 


3.  FLOWS  WITH  VORTICITY 

In  order  to  illustrate  the  effects  of  vorticity  two  sets  of  calculations 
(a  =  0.4,  8=1  end  a  =  0.5,  8  =  0.75)  were  investigated  at  inlet  Mach  numbers 
of  0.2  and  0.3.  "he  results  of  the  flow  field  calculations  are  presented  in 
Figs.  22  to  27.  The  distortion  of  axial  velocity  profile  due  to  the  preserce 
of  the  non-free  vortex  swirl  velocities  is  very  evident  in  the  figures.  Note, 
for  example  the  distorted  outlet  profiles  far  downstream.  Thus,  for  flows  with 
seemingly  low  approach  Mach  numbers,  the  distortion  of  the  profile  can  lead  to 
large  meridional  Mach  numbers  occurring  in  the  field.  The  very  strong  effect 
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figure  16.  Tangential  Velocity  Profiles  Immediately 
Following  Stator;  II  =  0.2,0. 3;  a  =  0, 


Figure  17.  Axial  Velocity  Profiles  at  Four  Axial  Stations; 


Figure  18.  Axial  Velocity  Profiles  at  Four  Axial  Stations; 


Figure  19.  Tangential  Velocity  Profiles  Immediately 
Following  Stator;  M  =  0.4, 0.5;  a  =  0,  £ 


Figure  20.  Axial  Velocity  Profiles  at  Four  Axial  Stations; 


Figure  21.  Axial  Velocity  Profiles  at  Four  Axial  Stations; 


Figure  22.  Tangential  Velocity  Profiles  Imnediately 
Following  Stator;  M  =  0.2,0. 3;  a  -  0.4, 
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Figure  25.  Tangential  Velocity  Profiles  Immediately  Following 
Stator;  M  =  0.2, 0.3;  a  «  0.5,  3  =  0.75 


Figure  27.  Axial  Velocitv  Profiles  at  Four  Axial  Stations; 


of  vorticity  can  be  noted  by  observing  the  values  of  the  axial  velocity  at  the 
hub  at  downstream  infinity  for  the  two  cases  at  M0  =  0.2.  Thus  when  a  ■  0.4, 
Wh(»)  =  1.67,  whereas  when  a  =  0.5,  W^(°°)  =  1.85.  The  latter  case  hence  leads 

more  rapidly  to  the  approach  to  supersonic  Mach  numbers.  This  is  true  even 
though  the  free  vortex  component  of  swirl  is  reduced  in  the  a  =  0.5  case  suffi¬ 
ciently  to  cause  the  average  outlet  enthalpies  to  be  equal. 

The  dotted  portions  of  the  curves  appearing  on  Figs.  24  and  27  indicate 
regions  where  the  computer  program  indicated  the  failure  of  the  density  itera¬ 
tion  program.  In  such  areas  the  density  is  estimated  by  the  calculational  pro¬ 
cedure  outlined  in  SECTION  IV-4.  Generally  speaking  the  number  of  mesh  points 
indicating  that  this  procedure  has  been  invoked  decreases  with  successive 
iterations.  In  these  cases  the  number  of  iterations  was  limited  to  ten,  so  it 
is  quite  possible  that  the  few  mesh  points  exhibiting  this  behavior  could  be 
removed  with  further  iterations.  As  it  is,  the  Meridional  Mach  number  exis  ing 
at  the  hub  at  downstream  infinity,  calculated  utilizing  the  density  ratio  ob¬ 
tained  by  the  alternate  calculation  procedure  is  .6^4  in  Fig  24  and  .975  in 
Fig.  27.  It  should  be  noted  that  the  values  chosen  for  a  (0.4  and  0.5)  repre¬ 
sent  very  high  levels  of  vorticity  (with  consequent  large  variations  in  stagna¬ 
tion  enthalpy)  compared  to  traditional  turbomachine  practice,  and  as  such  pro¬ 
vide  a  demanding  test  for  the  computer  program. 

4.  THE  EFFECT  OF  MESH  SIZE 


In  order  to  investigate  the  effect  of  Mesh  size,  some  very  highly  swirling 
free  vortex  flows  (8  =  6.66)  in  an  annulus  of  modest  tip  to  hub  ratio  (R  =  10/7) 
were  considered.  These  values  of  6  and  R  were  taken  in  order,  also,  to  compare 
the  results  to  those  of  Hawthorne  and  Ringrose(10)  (See  also  SECTION  VII).  The 
trends  observed  are  very  similar  to  those  reported  in  Ref.  10,  the  only  real 
differences  arising  because  of  slight  differences  in  the  value  of  8  assumed 
(6.66  versus  approximately  6.1  for  Ref.  10),  and  from  the  much  different  value 
of  wheel  speed  taken  here  (ft  =  1  versus  ft  =  2.61  for  Ref.  10). 

The  geometry  considered  has  parallel  walls,  with  the  stator  to  tu'or 
spacing  being  0.14.  (Thus  approximately  one  third  of  the  blade  height  R-l) . 

The  mesh  was  constructed  to  have  even  x  spacing  between  the  stator  and  rotor. 
Upstream  from  the  stator  and  dow-stream  from  the  rotor  fcne  x  spacing  of  the 
mesh  points  increases  in  the  ratio  1.25  for  each  succeeding  mesh  point.  In 
both  coarse  and  fine  meshes  there  are  eight  streamline  mesh  points.  In  the 
fine  mesh,  there  are  twelve  x  mesh  points  between  stator  and  rotor,  with  fifteen 
further  mesh  points  both  upstream  of  the  stator  and  downstream  of  the  rotor  for 
a  total  of  forty  two  x  mesh  points.  The  coarse  mesh  has  eight  mesh  points 
between  stator  and  rotor  with  eight  further  mesh  points  both  upstream  of  the 
stator  and  downstream  of  the  rotor  for  a  total  of  twenty  four  x  mesh  points. 

When  the  stator  alone  is  considered  the  mesh  point  spacing  is  the  same  as 
that  described  above.  Figs.  28-31  show  the  behavior  of  the  axial  velocity  on 
hub  and  tip  throughout  the  flow  field  length.  The  location  of  the  mesh  points 
is  indicated  in  the  figures. 

It  is  interesting  to  note,  again,  the  large  effect  of  'ompressibility  for 
cases  with  such  high  swirl  loading.  Thus,  though  the  entrance  Mach  number  is 
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Figure  28.  Coarse  Mesh  Results  for  Stator  Only 


Its,  Stator  and  Rotor 


Figure  31.  Fine  Mesh  Results,  Stator  and  Rotor 


only  Mq  =  0.182,  the  meridional  Mach  number  at  the  hub  just  following  the  stator 
is  0.497.  The  total  flow  Mach  number  at  this  location  is  1.54. 

It  can  be  seen  by  comparison  of  the  figures,  that  the  effect  of  decreased 
mesh  size  is  very  small  here.  (Careful  perusal  of  the  tabular  output  indicates 
changes  in  the  predicted  velocities  at  comparable  locations  of  only  approxi¬ 
mately  one  percent).  This  suggests  that  for  most  calculational  purposes  eight 
streamwise  mesh  points  and  twenty  four  x-wise  mesh  points  would  produce  suffi¬ 
cient  computational  accuracy. 

5.  THE  EFFECT  OF  FINITE  ASPECT  RATIO 


As  a  first  approximation  to  the  effect  of  finite  aspect  ratio  the  flow  of 
fluid  through  a  stator  only  was  considered.  The  first  calculation  (Fig.  32) 
considered  a  single  disc,  whereas  the  second  calculation  (Fig.  33)  considered  the 
blade  to  be  replaced  by  two  actuato1*  discs  such  that  62.5%  of  the  loading  was  on 
the  first  disc  and  37.5%  of  the  loading  on  the  second.  Figs.  30  and  31  show 
the  resultant  axial  velocity  profiles  existing  at  the  position  just  preceding 
the  upstream  disc  and  just  following  the  downstream  disc.  For  the  case  of  the 
single  disc,  the  x  =  0  origin  was  located  so  that  the  single  disc  was  at  the 
appropriate  center  of  moment. 

Comparison  of  the  profiles,  particularly  at  the  downstream  location,  indi¬ 
cate  as  expected,  that  the  flow  following  the  single  disc  has  proceeded  more 
towards  its  asymptotic  form  than  has  that  which  has  just  departed  the  second  of 
the  two  discs.  It  is  obvious  that  this  process  may  be  continued  to  a  larger 
number  of  actuator  discs  if  a  more  accurate  description  of  the  effects  of  dis¬ 
tributed  blade  forces  is  desired. 
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Figure  32.  Single  Disc  Approximation  to  Stator 


SECTION  VII 


COMPARISON  WITH  ALTERNATIVE  CALCULATIONAL  TECHNIQUES 


Many  calculational  techniques  for  throughflow  theory  are  available  in  the 
literature  today,  but  we  will  here  limit  our  attention  to  those  examples  that 
consider  compressible  flows. (7_1°) 

Creveling  and  CarmodyO)  provide  a  program  written  to  obtain  the  desired 
output  variables  of  a  turbomachine  row  or  rows  in  terms  of  input  design  para¬ 
meters.  The  study  is  not  really  directly  relatable  to  the  present  study,  be¬ 
cause  the  detailed  flow  information  involved  in  describing  the  flow  approach  to, 
and  departure  from,  a  blade  row  is  not  obtained.  Rather,  a  "quasi  radial 
equilibrium"  solution  is  obtained  for  the  conditions  between  blade  rows.  In 
this  quasi  radial  equilibrium  solution  the  effect  of  streamline  curvature  at  the 
station  investigated  is  incorporated  directly  by  taking  the  supplied  flow  path 
geometry  (Mod  I)  or  by  establishing  the  flow  path  geometry  by  computation 
(Mod  II).  It  is  pointed  out  by  the  authors  that  irregular  geometry  can  be  pro¬ 
duced  by  Mod  II  computation,  resulting  in  severe  streamline  curvature  effects. 

In  common  with  all  the  example  techniques  studied  (and  in  common  with  this 
report)  the  computation  is  limited  to  flows  in  which  the  meridional  Mach  number 
is  less  than  unity.  By  simplifynng  the  fluid  mechanical  description,  the 
authors  have  made  it  possible  to  analyze  the  machine  behavior  in  terms  of  de¬ 
sign  input  variables  (total  pressure  profile  at  each  rotor  exit,  and  limiting 
values  of  Rotor  tip  diffusion  factor,  Stator  hub  diffusion  factor,  Stator  hub 
Mach  number,  Rotor  hub  relative  exit  angle  and  Rotor  tip  exit  whirl  velocity). 
The  output  of  the  program  is  the  desired  flow  field  information.  (Velocities, 
temperatures  flow  angles,  etc.) 

The  study  by  Marsh(8)  is  directed  to  the  solution  of  the  "direct  problem" 
of  throughflow  theory  wherein  the  blade  geometry  is  prescribed  and  the  resulting 
flow  field  is  to  be  determined.  As  a  result  actuator  disc  theory  is  not  used, 
but  rather  the  "average  forces"  of  the  blades  are  related  to  the  rate  of 
change  of  angular  momentum  of  the  fluid  and  the  blade  surface  geometry.  The 
effect  of  blade  thickness  is  included  by  adding  a  factor,  B,  in  the  continuity 
equation  that  is  simply  the  ratio  of  the  circumferential  width  of  the  blade 
passage  to  the  blade  pitch.  The  effects  of  imperfect  flow  can  be  included 
through  the  generation  of  entropy  term,  though  the  related  drag  effect  upon  the 
equation  of  motion  is  not  included. 

It  is  interesting  to  note  that  in  common  with  Wu,(1J)  Marsh  elected  to 
jiJve  and  store  the  equation  for  density  (our  Eq  (13))  rather  than  iterate  for 
the  solution  at  each  point.  In  the  examples  considered  in  the  present  report, 
after  the  first  pass  for  streamline  location  the  number  of  iterations  required 
for  six  figure  density  accuracy  were  usually  only  one  or  two.  The  authors  of 
this  report  feel  that  when  such  accuracy  is  desired  it  is  more  efficient  to  use 
the  iterative  technique  than  to  store  the  information. 

The  actual  computational  process  used  by  Marsh  involves  solving  a  quasi- 
linear  equation  for  T  in  terms  of  an  assumed  known  density  field.  After  solu¬ 
tion  of  the  quasi-linear  equation  the  density  estimate  is  updated  and  the  pro¬ 
cess  repeated  till  convergence.  Finite  difference  methods  are  used  throughout. 
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The  examples  given  by  Marsh  correspond  to  very  low  Mach  number  flows 
(speeds  of  less  than  100  fps) ,  so  it  is  difficult  to  determine  with  certainty 
the  reliability  of  the  convergence  of  the  process  when  applied  to  the  analysis 
of  flows  in  which  large  meridional  Mach  numbers  exist.  Marsh  mentions  that  an 

V  -  T, 

—  3  —  tf  k. 

accuracy  of  10  to  10  in  the  streamline  position  estimate,  — ^ -  is  normal- 

K 

ly  obtained  within  10  to  20  iterations. 

Schroder  and  Schuster(9)  employ  a  finite  difference  technique  to  solve  the 
compressible  throughflow  problem  in  a  parallel  walled  annulus  containing  one  or 
two  actuator  discs.  Their  formulation  is  "quasi-direct"  in  the  sense  that  they 
prescribe  the  jump  in  tangential  velocity  (and  hence  in  enthalpy  when  the  disc 
is  a  rotor)  in  terms  of  the  spacial  (r)  position  rather  than  in  terms  of  the 
value  of  the  streamfunction.  By  so  doing,  they  are  forced  to  resort  to  itera¬ 
tive  corrections  to  the  equations  describing  the  flow  downstream  of  a  rotor 
because  the  local  value  of  stagnation  enthalpy  at  a  given  point  (r,z)  must  be 
updated  as  the  value  of  the  streamfunction  at  the  disc  and  at  the  point  (r,z), 
change.  As  with  Ref.  8  and  the  present  report,  the  density  is  calculated  by 
iteratively  updating  the  streamline  position  estimate. 

In  common  with  the  present  report,  Schroder  and  Schuster  consider  only 
radial  blades.  The  main  thrust  of  the  paper  of  Ref.  9  seems  to  be  in  the  direc¬ 
tion  of  investigating  the  jurne  conditions  at  the  actuator  disc.  It  appears 
that  their  program  is  the  most  developed  with  regard  to  capability  of  handling 
high  meridional  Mach  numbers  of  the  papers  reviewee  here,  but  the  parallel 
walled  geometry  is  unfortunately  restrictive.  With  regards  to  convergence, 
Schroder  and  Schuster  state  that  "as  a  rule,  35  to  70  complete  iteration  cycles 
will  be  necessary  to  arrive  at  convergence,  each  cycle  passing  over  all  nodes." 

The  final  paper  considered  here,  that  by  Hawthorne  and  Ringrose(10)  is  the 
most  restrictive  physically  in  that  it  considers  only  free  vortex  flow  through 
actuator  discs  contained  in  a  parallel  walled  annulus.  The  solutions,  however, 
are  obtained  in  analytical  form  by  first  linearizing  the  equations.  It  is 
apparent  that  the  linearizing  approximation  requires  the  presence  of  "small" 
meridional  Mach  numbers  in  order  to  give  valid  results.  The  paper,  however,  is 
very  useful  for  several  reasons.  Most  importantly  the  very  rapid  calculation 
of  desired  examples  allows  the  easy  physical  interpretation  of  the  density  ef¬ 
fects.  In  addition,  by  considering  a  free  vortex  flow  field,  no  velocity 
perturbations  arise  from  rotational  effects  so  the  effects  of  density  varia¬ 
tions  are  effectively  isolated.  Fina’ly,  the  restriction  to  low  meridional 
Mach  numbers  is  not  as  restrictive  as  at  first  might  seem  apparent,  because 
flows  with  very  large  swirl  Mach  numbers  (>  1)  can  be  considered.  Such  flows 
exhibit  large  compressibility  effects  but  do  not  necessarily  contain  large 
meridional  Mach  numbers.  Several  computational  examples  have  been  given  already 
in  SECTION  VI  in  which  results  of  the  present  report  calculations  were  compared 
to  the  results  of  Ref.  10. 
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SECTION  VIII 


SUMMARY  AND  CONCLUSIONS 


This  report  has  given  the  statement  of  a  new  variational  formulation  of  the 
compressible  throughflow  problem.  A  detailed  development  of  a  finite  element 
analysis  and  presentation  of  an  associated  computer  program  is  included,  together 
with  many  example  calculations. 

It  is  important  to  note  that  the  variational  principle  is  valid  for  a  very 
general  flow  field.  Thus,  compressibility,  entropy  variations  and  real  gas  ef¬ 
fects  may  all  be  included  in  the  analysis.  It  is  assumed  that  viscosity  is  zero 
throughout  the  flow  field,  an  assumption  that  implies  that  any  entropy  variations 
arise  within  the  blade  rows.  This  latter  assumption  is,  of  course,  consistent 
with  conventional  throughflow  theory. 

The  analytical  and  numerical  forms  presented  herein  have  been  restricted  to 
isentropic  flows  and  to  flows  of  a  calorically  perfect  gas.  This  latter  assump¬ 
tion  appears  primarily  in  the  subsidiary  expression  for  the  density.  It  is 
emphasized  again,  however,  that  further  development  of  the  concept  could  lead  to 
inclusion  of  both  entropy  variations  and  real  gas  effects. 

The  authors  feel  that  the  utilization  of  finite  element  techniques  in  this 
context  promises  to  lead  to  the  efficient  calculation  of  very  complex  flow  fields 
and  hope  that  this  report  has  made  a  significant  step  in  this  direction. 
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